Over-extension of the fault models

[Geo-Explore]

Hello! I have been using Gempy and it had worked great so far!

I have successfully constructed a series of fault models, but a problem has arisen at the moment, due to the large extent of my study area and the fact that different faults are located in different areas, I have set up a large range of models to ensure that my orientation data and observations are included as much as possible.

However, this has led to a problem where the modeled faults have extended themselves to a certain extent, and faults have intersected that should not have intersected in the first place. I'm wondering if there is a way to solve this problem - how to construct multiple faults at large regional scales and have them not over-extend, and only construct the parts that have data.

I'd appreciate any help/hints what path to follow next! Thanks in advance!

[AlexanderJuestel]

See #675 and #787 for example. It is currently not fully implemented but there may be ways to trick the system.

Let us know if you have any other issues!

[Geo-Explore]

@AlexanderJuestel Great! Thank you for your response!

[AlexanderJuestel]

I will also think of a little workaround one of these days that I would implement in our GemGIS package

[Geo-Explore]

@AlexanderJuestel That's great! Looking forward to your updates!

[AlexanderJuestel]

The first idea (for limiting the faults only to the interface points): Creating a plane perpendicular to the fault orientation at the last interface. Then clipping the fault by the plane. The second plane illustrates a plane that was translated parallel to the strike of the fault and could clip the fault also at that point.

These tasks would be performed in post-processing.

[Geo-Explore]

@AlexanderJuestel I apologize for the late reply. That sounds great. But I have a question, can this plane used for truncated faults be subsequently deleted? Because I only want to keep the fault plane in the geologic model and don't want to include other extraneous planes.

[AlexanderJuestel]

Yes! I am currently creating a function that will only return the clipped fault. Showing the planes was only for illustration now ;)

I should have a proper solution by tomorrow!

[Geo-Explore]

That's cool！

---- Replied Message ---- | From | @.> | | Date | 08/17/2023 21:16 | | To | @.> | | Cc | @.>@.> | | Subject | Re: [cgre-aachen/gempy] Over-extension of the fault models (Issue #821) |

Yes! I am currently creating a function that will only return the clipped fault. Showing the planes was only for illustration now ;)

I should have a proper solution by tomorrow!

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: @.***>

[AlexanderJuestel]

@Geo-LiuDaLi,

I think I got something pretty solid. The gray fault is a fault modeled by GemPy which is extended beyond the interface points. The red dots are the interface points of the faults (first and last interface point). The reddish part is the clipped fault. On the right side at the last interface and on the left side with a buffer of 250 m along the strike of the fault.

What needs to be fixed first is https://github.com/cgre-aachen/gemgis/issues/288 (fix will be merged for GemGIS 1.1) and I requested a little feature in https://github.com/cgre-aachen/gempy/pull/822 that I would currently need to achieve result.

This would be my code but as mentioned, there are two fixes needed to make it running (it runs locally)

from typing import Union, List
import geopandas as gpd
# from gemgis.visualization import create_depth_maps_from_gempy

def clip_fault_of_gempy_model(geo_model,
                              fault: str,
                              which: str = 'first',
                              buffer_first: Union[int, float] = None,
                              buffer_last: Union[int, float] = None,
                              i_size: Union[int, float] = 1000,
                              j_size: Union[int, float] = 1000,
                              invert_first: bool = True,
                              invert_last: bool = False) -> Union[pv.core.pointset.PolyData, List[pv.core.pointset.PolyData]]:
    """
    Clip fault of a GemPy model.

    Parameters
    __________
        geo_model : gp.core.model.Project
            GemPy Model containing the faults.
        fault : str
            String or list of strings containing the name of faults to be clipped, e.g. ``faults='Fault1'``-
        which : str, default: ``'first'``
            Parameter to decide which end of the faults to clip. Options include ``'first'``, ``'last'``, or both, e.g. ``'which='first'``.
        buffer_first : Union[int, float]
            Int or float value or list of values to clip the fault/s behind the first interface point, e.g. ``'buffer_first=500'``.
        buffer_last : Union[int, float]
            Int or float value or list of values to clip the fault/s behind the last interface point, e.g. ``'buffer_last=500'``.
        i_size: Union[int, float]
            Size of the plane in the i direction.
        j_size: Union[int, float]
            Size of the plane in the j direction.
        invert_first : bool, default: ``'True'``
            Invert clipping for first plane.
        invert_last : bool, default: ``'False'``
            Invert clipping for second plane.

    Returns
    _______
        pv.core.pointset.PolyData
            Clipped faults

    .. versionadded :: 1.1

    See also
    ________
        create_plane_from_interface_and_orientation : Create PyVista plane from GemPy interface and orientations DataFrames.
        translate_clipping_plane : Translate clipping plane.

    Example
    _______

    """
    # Checking that the fault is provided as string 
    if not isinstance(fault, str):
        raise TypeError('Faults must be provided as one string for one fault ')

    # Checking that the fault is a fault of the geo_model
    if isinstance(fault, str):
        if not fault in geo_model.surfaces.df['surface'][geo_model.surfaces.df['isFault']==True].tolist():
            raise ValueError('Fault is not part of the GemPy geo_model')

        # Getting the fault DataFrames
        fault_df_interfaces = gp.get_data(geo_model, 'surface_points')[gp.get_data(geo_model, 'surface_points')['surface']==fault].reset_index(drop=True)
        fault_df_orientations = gp.get_data(geo_model, 'orientations', numeric=True)[gp.get_data(geo_model, 'orientations', numeric=True)['surface']==fault].reset_index(drop=True)

    # Checking that the parameeter which is of type string or list of strings
    if not isinstance(which, str):
        raise TypeError('The parameter "which" must be provided as string. Options for each fault include "first", "last", or "both"')

    # Checking that the correct values are provided for the parameter which
    if isinstance(which, str):
        if not which in ['first', 'last', 'both']:
            raise ValueError('The options for the parameter "which" include "first", "last", or "both"')

    # Checking that the i size is of type int or float
    if not isinstance(i_size, (int, float)):
        raise TypeError('i_size must be provided as int or float')

    # Checking that the j size is of type int or float
    if not isinstance(j_size, (int, float)):
        raise TypeError('j_size must be provided as int or float')

    # Extracting depth map
    mesh = create_depth_maps_from_gempy(geo_model, 
                                        surfaces=fault)

    # Getting the first interface points
    if which == 'first':

        fault_df_interfaces_selected = fault_df_interfaces.iloc[0:1].reset_index(drop=True)

        # Creating plane from DataFrames
        plane, azimuth = create_plane_from_interface_and_orientation_dfs(df_interface=fault_df_interfaces_selected, 
                                                                         df_orientations=fault_df_orientations,
                                                                         i_size=i_size,
                                                                         j_size=j_size)
        # Translating Clipping Plane
        if buffer_first:

            # Checking that buffer_first is of type int or float
            if not isinstance(buffer_first, (int, float)):
                raise TypeError('buffer_first must be provided as int or float')

            plane = translate_clipping_plane(plane=plane, 
                                             azimuth=azimuth, 
                                             buffer=buffer_first)
        # Clipping mesh
        mesh[fault][0] = mesh[fault][0].clip_surface(plane, 
                                                     invert=invert_first)

    # Getting the last interface points
    elif which == 'last':

        fault_df_interfaces_selected = fault_df_interfaces.iloc[-1:].reset_index(drop=True)

        # Creating plane from DataFrames
        plane, azimuth = create_plane_from_interface_and_orientation_dfs(df_interface=fault_df_interfaces_selected, 
                                                                         df_orientations=fault_df_orientations,
                                                                         i_size=i_size,
                                                                         j_size=j_size)

        # Translating Clipping Plane
        if buffer_last:

            # Checking that buffer_last is of type int or float
            if not isinstance(buffer_last, (int, float)):
                raise TypeError('buffer_last must be provided as int or float')

            plane = translate_clipping_plane(plane=plane, 
                                             azimuth=azimuth, 
                                             buffer=buffer_last)

        # Clipping mesh
        mesh[fault][0] = mesh[fault][0].clip_surface(plane, 
                                                     invert_last)

    if which == 'both':

        # First point
        fault_df_interfaces_selected = fault_df_interfaces.iloc[0:1].reset_index(drop=True)

        # Creating plane from DataFrames
        plane1, azimuth1 = create_plane_from_interface_and_orientation_dfs(df_interface=fault_df_interfaces_selected, 
                                                                           df_orientations=fault_df_orientations,
                                                                           i_size=i_size,
                                                                           j_size=j_size)
        # Translating Clipping Plane
        if buffer_first:
            plane1 = translate_clipping_plane(plane=plane1, 
                                              azimuth=azimuth1, 
                                              buffer=buffer_first)

        # Last Point
        fault_df_interfaces_selected = fault_df_interfaces.iloc[-1:].reset_index(drop=True)

        # Creating plane from DataFrames
        plane2, azimuth2 = create_plane_from_interface_and_orientation_dfs(df_interface=fault_df_interfaces_selected, 
                                                                           df_orientations=fault_df_orientations,
                                                                           i_size=i_size,
                                                                           j_size=j_size)

        # Translating Clipping Plane
        if buffer_last:
            plane2 = translate_clipping_plane(plane=plane2, 
                                              azimuth=azimuth2, 
                                              buffer=-buffer_last)

        # Clipping mesh
        mesh[fault][0] = mesh[fault][0].clip_surface(plane1, 
                                                     invert=invert_first).clip_surface(plane2, 
                                                                                       invert=invert_last)

    return mesh

def create_plane_from_interface_and_orientation_dfs(df_interface: pd.DataFrame,
                                                    df_orientations: pd.DataFrame, 
                                                    i_size: Union[int, float] = 1000,
                                                    j_size: Union[int, float] = 1000) -> pv.core.pointset.PolyData:
    """ 
    Create PyVista plane from GemPy interface and orientations DataFrames.

    Parameters
    __________
        df_interface : pd.DataFrame
            GemPy Pandas DataFrame containing the interface point for the plane creation.
        df_orientations : pd.DataFrame
            GemPy Pandas Dataframe containing the orientations for the plane creation.
        i_size: Union[int, float]
            Size of the plane in the i direction.
        j_size: Union[int, float]
            Size of the plane in the j direction.

    Returns
    _______
        plane : pv.core.pointset.PolyData
            Plane for clipping the fault. 
        azimuth : Union[int, float]
            Azimuth of the fault. 

    .. versionadded:: 1.1

    See also
    ________
        clip_fault_of_gempy_model : Clip fault of a GemPy model.
        translate_clipping_plane : Translate clipping plane.

    Example
    _______

    """
    # Checking that the interface DataFrame is a DataFrame
    if not isinstance(df_interface, pd.DataFrame):
        raise TypeError('Interface must be provided as Pandas DataFrame')

    # Checking that the orientations DataFrame is a DataFrame
    if not isinstance(df_orientations, pd.DataFrame):
        raise TypeError('Orientations must be provided as Pandas DataFrame')

    # Checking that the i size is of type int or float
    if not isinstance(i_size, (int, float)):
        raise TypeError('i_size must be provided as int or float')

    # Checking that the j size is of type int or float
    if not isinstance(j_size, (int, float)):
        raise TypeError('j_size must be provided as int or float')

    # Creating GeoDataFrame from interface
    gdf_interface = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=df_interface['X'],
                                                                 y=df_interface['Y']), 
                                     data=df_interface)

    # Creating GeoDataFrame from orientations
    gdf_orientations = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=df_orientations['X'],
                                                                    y=df_orientations['Y']), 
                                        data=df_orientations)

    # Finding nearest orientation to the respective interface to set the orientation of the plane
    gdf_orientations_nearest = gpd.sjoin_nearest(gdf_interface, 
                                                 gdf_orientations)

    # Extracting azimuth for clipping plane
    azimuth = gdf_orientations_nearest['azimuth'][0]

    # Extracting center of clipping plane
    center = df_interface[['X', 'Y', 'Z']].values[0]

    # Creating clipping plane, direction is created from the orientation of the fault.
    plane = pv.Plane(center=center, 
                     direction=(np.cos(np.radians(azimuth)), 
                                np.sin(np.radians(azimuth)), 
                                0.0), 
                     i_size=i_size, 
                     j_size=j_size)

    return plane, azimuth

def translate_clipping_plane(plane: pv.core.pointset.PolyData,
                             azimuth: Union[int, float, np.int64],
                             buffer: Union[int, float]) -> pv.core.pointset.PolyData:
    """
    Translate clipping plane.

    Parameters
    __________
        plane : pv.core.pointset.PolyData
            Clipping Plane.
        azimuth : Union[int, float, np.int64]
            Orientation of the Fault.
        buffer : Union[int, float, type(None)]
            Buffer to translate the clipping plane along the strike of the fault.

    Returns
    _______
        pv.core.pointset.PolyData
            Translated clipping plane.

    .. versionadded:: 1.1

    See also
    ________
        create_plane_from_interface_and_orientation : Create PyVista plane from GemPy interface and orientations DataFrames.
        clip_fault_of_gempy_model : Clip fault of a GemPy model.

    Example
    _______

    """
    # Checking that the plane is of type PyVista PolyData
    if not isinstance(plane, pv.core.pointset.PolyData):
        raise TypeError('The clipping plane must be provided as PyVista PolyData')

    # Checking that the azimuth is of type int or float
    if not isinstance(azimuth, (int, float, np.int64)):
        raise TypeError('The azimuth must be provided as int or float')

    # Checking that the buffer is of type int or float
    if not isinstance(buffer, (int, float, type(None))):
        raise TypeError('The buffer must be provided as int or float')

    # Calculating translation factor in X and Y Directio
    x_translation = -np.cos(np.radians(azimuth))*buffer
    y_translation = -np.sin(np.radians(azimuth))*buffer

    # Translating plane
    plane = plane.translate((x_translation*np.cos(np.radians(azimuth)), 
                                 y_translation*np.sin(np.radians(azimuth)), 
                                 0.0))

    return plane

[Japhiolite]

short question: This only affects the visual presentation of the fault plane, right? Not the distortion of the scalar field, i.e. the displacement caused by the fault? We could think of directly implementing this in gempy instead of gemgis and use it if a user sets the fault to finite.

[AlexanderJuestel]

@Japhiolite, that is correct, yes. I am basically just clipping the fault planes.

[AlexanderJuestel]

@Geo-LiuDaLi Have a look at this notebook here and let me know if the solution works for you. You will have to install GemGIS Version 1.1 using pip install gemgis==1.1 for that https://gemgis.readthedocs.io/en/latest/getting_started/tutorial/68_Creating_Finite_Faults_with_GemGIS.html

[Geo-Explore]

@AlexanderJuestel Thank you very much! I'll get back to you as soon as I can.

[Geo-Explore]

Hi @AlexanderJuestel ！ I tried to use this solution to build finite faults in my model, after I added the appeal fix, I can build the base fault model properly, but a new problem appeared and does not enable clipping of faults, can you help me? Attached is the error message, please feel free to contact me if you need any other information! Thanks in advance!

[AlexanderJuestel]

@Geo-LiuDaLi I thought I had fixed that issue in GemGIS 1.1. Which versions of PyVista and NumPy are you using?

[Geo-Explore]

Hi@AlexanderJuestel I'm using: NumPy : 1.21.6 PyVista : 0.34.2 GemGIS : 1.1.0 GemPy : 2.2.12

[AlexanderJuestel]

Please try upgrading GemPy to version 2.3 (pip install gempy==2.3) and PyVista (pip install pyvista==0.41.1) and NumPy (pip install numpy==1.25.2)

[Geo-Explore]

Okay. Thank you very much. I'll try an updated version as soon as I can.

[AlexanderJuestel]

@Geo-LiuDaLi Have you had the chance to test the clipping of faults?

[Geo-Explore]

Hi @AlexanderJuestel ~ I've tried to test the clipping feature for faults, but it's rendering may not be quite what I expected. Upon testing I've found that it currently only seems to change the presentation of existing models visually, and doesn't really enable clipping of faults. I'm not quite sure if I'm doing something wrong somewhere to the extent that it's causing this problem. Thanks for the help anyway! :）

[AlexanderJuestel]

Yes, you are correct! It is just changing the appearance of the fault. It does not interfere with the GemPy modeling itself

[Geo-Explore]

Alright, it doesn't seem to be quite what I was expecting. I wanted to achieve a clipping of the fracture model itself and not just change the appearance of the fracture, in other words, I wanted to build truly finite faults. Anyway, thank you very much for your help, and if you have any new ideas about the construction of finite faults please feel free to contact me, thanks!:)

[Leguark]

Finite fault are not officially supported yet. There are some tests of the first implementations for brave people but nothing exposed on the API