draw a geological faults (break down in surface of rock formations

I do a trail of following MATLAB code : to clasfying the geological faults which are broken on rock formations represnting by inverted depths ( z_inv_Bouguer) whch can reading from the annexed Excel file : I need to draw faults like in annexed figure dependenig on basis of angles (Theta) drived from gradient of each curve.

N = 9; % Number of stacked layers

M = 201; % Number of digitized points

xc = data(:, 1); % x-axis direction

z_inv_Bouguer = data(:,3:end) ; % z_inv_Bouguer curves

% Initialize arrays to store fault locations based on angle classifications

normal_fault_locations = [];

reverse_fault_locations = [];

strike_slip_fault_locations = [];

other_fault_locations = [];

for i = 1:N

% Calculate gradient of z_inv_Bouguer curves

for j = 1:M

Theta = rad2deg(atan(z_inv_Bouguer(j,i) ./ xc(j)));

% Define thresholds for classifying faults

normal_fault_threshold = 90; % Adjust as needed

reverse_fault_threshold = 45; % Adjust as needed

% Classify faults based on angle

if abs(Theta) <= normal_fault_threshold

normal_fault_locations = [normal_fault_locations; [xc(j),...

z(i)]];

elseif abs(Theta) > reverse_fault_threshold

reverse_fault_locations = [reverse_fault_locations; [xc(j),...

z(i)]];

elseif abs(Theta) > normal_fault_threshold...

&& abs(Theta) <= reverse_fault_threshold

strike_slip_fault_locations = [strike_slip_fault_locations;

[xc(j), z(i)]];

else

other_fault_locations = [other_fault_locations; [xc(j), z(i)]];

end

end

end

% Plotting results

figure;

% Plot z_inv_Bouguer curves

for i = 1:N

plot(xc, z_inv_Bouguer(:, i), 'LineWidth', 1);

hold on;

end

% Plot normal fault locations

if ~isempty(normal_fault_locations)

scatter(normal_fault_locations(:, 1),...

normal_fault_locations(:, 2),...

'r', 'filled', 's', 'DisplayName', 'Normal Faults');

hold on;

end

% Plot reverse fault locations

if ~isempty(reverse_fault_locations)

scatter(reverse_fault_locations(:, 1)...

, reverse_fault_locations(:, 2)...

, 'b', 'filled', 'o', 'DisplayName', 'Reverse Faults');

hold on;

end

% Plot strike-slip fault locations

if ~isempty(strike_slip_fault_locations)

scatter(strike_slip_fault_locations(:, 1)...

, strike_slip_fault_locations(:, 2)...

, 'g', '+', 'DisplayName', 'Strike-Slip Faults');

hold on;

end

% Plot other fault locations (if any)

if ~isempty(other_fault_locations)

scatter(other_fault_locations(:, 1)...

, repmat(min(z), size(other_fault_locations, 1), 1)...

, 'm', 'filled', 'v', 'DisplayName', 'Other Faults');

hold on;

end

set(gca, 'YDir', 'reverse');

title(' z_inv_Bouguer with Classified Faults');

xlabel('Profile Points-xc');

ylabel('Depth (km)');

legend('show');

grid on;

%Thanks in advance

Answer 1

[Suraj Kumar]

Hi Moustafa,

To classify and plot geographical faults using the gradient of inverted depth curves from an Excel file, follow these steps and the attached code snippets:

1. Load the data from the excel file, define the parameters and extract the values representing the inverted depth curves.

data = xlsread('z_inv_Bouguer.xlsx');  
N = 8
M = 201;  
xc = data(:, 1);   
z_inv_Bouguer = data(:, 3:end);  

2. Initialize arrays that will store the coordinates of the classified fault locations.

normal_fault_locations = [];
reverse_fault_locations = [];
strike_slip_fault_locations = [];
other_fault_locations = [];

3. Calculate the gradient of the “z_inv_Bouguer”curves by iterating through each layer and point, using differences in “z_inv_Bouguer”and “xc” values.

for i = 1:N
    for j = 2:M-1 
        dz = z_inv_Bouguer(j+1, i) - z_inv_Bouguer(j-1, i);
        dx = xc(j+1) - xc(j-1);
        Theta = rad2deg(atan(dz / dx));
        normal_fault_threshold = 45; 
        reverse_fault_threshold = -45; 

4. Classify the faults based on the angle theta and store them in the respective array.

if Theta > normal_fault_threshold
    normal_fault_locations = [normal_fault_locations; [xc(j), z_inv_Bouguer(j, i)]];
 elseif Theta < reverse_fault_threshold
    reverse_fault_locations = [reverse_fault_locations; [xc(j), z_inv_Bouguer(j, i)]];
 elseif Theta <= normal_fault_threshold && Theta >= reverse_fault_threshold
    strike_slip_fault_locations = [strike_slip_fault_locations; [xc(j), z_inv_Bouguer(j, i)]];
 else
    other_fault_locations = [other_fault_locations; [xc(j), z_inv_Bouguer(j, i)]];
 end

5. Plot the depth curves and the fault locations to indicate different fault types.

% Plot z_inv_Bouguer curves
for i = 1:N
    plot(xc, z_inv_Bouguer(:, i), 'LineWidth', 1);
    hold on;
end

Refer to the output for better understanding:

To learn more about “atan” and “rad2deg” functions in MATLAB, kindly refer to the following documentations:

• https://www.mathworks.com/help/matlab/ref/atan.html
• https://www.mathworks.com/help/matlab/ref/rad2deg.html

I hope this helps.