I for got the files
Need to put symbols on curve
4 views (last 30 days)
Show older comments
Moustafa Abedel Fattah
on 26 Mar 2024
Commented: Mathieu NOE
on 27 Mar 2024
I write the following Matlab code (R2014b) for detecting the geological fault (break down) and the iclination angle and type of fault and I need to put the symboles of fault as in the secobd figuer of the annexed fike - the RasGhari BouguetProfile_AB.xlsx is anned
Thanks in advance
%===========================================================%
clc
close all;
clear;
workspace;
%===========================================================%
% Constants
G = 0.0067;
delt_rohj = 0.050;
% Load data from Excel file
% data = xlsread('RasGhari BouguetProfile_AB.xlsx');
data = xlsread('RasGhariSepBougtProfiles.xlsx');
xc = data(:, 1); % x-coordinate
gBt = data(:, 2); % Bouguer gravity anomaly
% Calculate gradient of the gravity anomaly profile
gradient_gBt = abs(gradient(gBt));
% Estimate the gradient threshold automatically
gradient_threshold = median(abs(gradient_gBt)) ; % Can adjust by multiplier as needed
% Detect fault locations
fault_locations = find(abs(gradient_gBt) > gradient_threshold);
% Initialize arrays to store fault angles and types
fault_angles = zeros(size(fault_locations));
fault_types = cell(size(fault_locations));
% Parameters
G = 6.67430e-3; % Gravitational constant
delt_rohj = 0.050; % Density contrast (example value, adjust as needed)
for i = 1:numel(fault_locations)
% Get indices for calculating fault properties
j = fault_locations(i);
if j == 1 || j == numel(xc)
continue; % Skip if fault is at the boundaries
end
% Calculate fault angle
C = abs(atan((gBt(j+1) - gBt(j-1)) / (xc(j+1) - xc(j-1)))); % Angle in Radians
B = rad2deg(C); % Angle in Degrees
% Calculate h (difference in gBt values)
h = gBt(i+1) - gBt(i);
g_fault_i = 2 * pi * G * delt_rohj * h *...
(1 + (1/pi) * atan(xc(i)./gBt(i) + cotd(B)) - (1/pi)...
* atan(xc(i)./gBt(i+1) + cotd(B)));
% Determine fault type
if B <= 90;
fault_type = 'Normal';
else
fault_type = 'Reverse';
end
% Store fault angle and type
fault_angles(i) = 90-B;
fault_types{i} = fault_type;
end
% Display results
disp('Detected Faults:');
disp('Location Angle (degrees) Type');
for i = 1:numel(fault_locations)
disp([num2str(xc(fault_locations(i))), ' ', num2str(fault_angles(i)), ' ', fault_types{i}]);
end
% Plotting
figure(1)
plot(xc, gBt, 'LineWidth', 1);
hold on;
grid on;
title('Bouguer Anomaly');
xlabel('Profile Points-xc');
ylabel('Gravity Anomaly');
legend('Bouguer Anomaly');
figure(2)
plot(xc, gBt, 'LineWidth', 1);
hold on;
grid on;
title('Bouguer Anomaly');
xlabel('Profile Points-xc');
ylabel('Gravity Anomaly');
% Add symbols for fault angles
for i = 1:numel(fault_locations)
x_fault = xc(fault_locations(i));
y_fault = gBt(fault_locations(i));
angle = fault_angles(i);
if strcmp(fault_types{i}, 'Normal')
symbol = 'o'; % Use circle for normal faults
else
symbol = 'x'; % Use cross for reverse faults
end
scatter(x_fault, y_fault, 100, 'Marker', symbol, 'MarkerEdgeColor', 'r', 'LineWidth', 2);
% text(x_fault, y_fault, sprintf('%0.1f°', angle), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
end
legend('Bouguer Anomaly', 'Faults', 'Location', 'best');
Accepted Answer
Mathieu NOE
on 26 Mar 2024
Like that ?
%===========================================================%
clc
close all;
clear;
workspace;
%===========================================================%
% Constants
G = 0.0067;
delt_rohj = 0.050;
% Load data from Excel file
% data = xlsread('RasGhari BouguetProfile_AB.xlsx');
data = xlsread('RasGhariSepBougtProfiles.xlsx');
xc = data(:, 1); % x-coordinate
gBt = data(:, 2); % Bouguer gravity anomaly
% Calculate gradient of the gravity anomaly profile
gradient_gBt = abs(gradient(gBt));
% Estimate the gradient threshold automatically
gradient_threshold = median(abs(gradient_gBt)) ; % Can adjust by multiplier as needed
% Detect fault locations
fault_locations = find(abs(gradient_gBt) > gradient_threshold);
% Initialize arrays to store fault angles and types
fault_angles = zeros(size(fault_locations));
fault_types = cell(size(fault_locations));
% Parameters
G = 6.67430e-3; % Gravitational constant
delt_rohj = 0.050; % Density contrast (example value, adjust as needed)
for i = 1:numel(fault_locations)
% Get indices for calculating fault properties
j = fault_locations(i);
if j == 1 || j == numel(xc)
continue; % Skip if fault is at the boundaries
end
% Calculate fault angle
C = abs(atan((gBt(j+1) - gBt(j-1)) / (xc(j+1) - xc(j-1)))); % Angle in Radians
B = rad2deg(C); % Angle in Degrees
% Calculate h (difference in gBt values)
h = gBt(i+1) - gBt(i);
g_fault_i = 2 * pi * G * delt_rohj * h *...
(1 + (1/pi) * atan(xc(i)./gBt(i) + cotd(B)) - (1/pi)...
* atan(xc(i)./gBt(i+1) + cotd(B)));
% Determine fault type
if B <= 90
fault_type = 'Normal';
else
fault_type = 'Reverse';
end
% Store fault angle and type
fault_angles(i) = 90-B;
fault_types{i} = fault_type;
end
% Display results
disp('Detected Faults:');
disp('Location Angle (degrees) Type');
for i = 1:numel(fault_locations)
disp([num2str(xc(fault_locations(i))), ' ', num2str(fault_angles(i)), ' ', fault_types{i}]);
end
% Plotting
figure(1)
plot(xc, gBt, 'LineWidth', 1);
hold on;
grid on;
title('Bouguer Anomaly');
xlabel('Profile Points-xc');
ylabel('Gravity Anomaly');
legend('Bouguer Anomaly');
figure(2)
plot(xc, gBt, 'LineWidth', 1);
hold on;
grid on;
title('Bouguer Anomaly');
xlabel('Profile Points-xc');
ylabel('Gravity Anomaly');
% Add symbols for fault angles
for i = 1:numel(fault_locations)
x_fault = xc(fault_locations(i));
y_fault = gBt(fault_locations(i));
angle = fault_angles(i);
if strcmp(fault_types{i}, 'Normal')
symbol = 'o'; % Use circle for normal faults
else
symbol = 'x'; % Use cross for reverse faults
end
scatter(x_fault, y_fault, 100, 'Marker', symbol, 'MarkerEdgeColor', 'r', 'LineWidth', 2);
% text(x_fault, y_fault, sprintf('%0.1f°', angle), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
end
% start and end indexes of fault locations
ind_end = find(diff(fault_locations)>1);
ind_start = [1; ind_end+1];
ind_all = [ind_start;ind_end];
for k = 1:numel(ind_all)
xx = xc(fault_locations(ind_all(k)));
yy = gBt(fault_locations(ind_all(k)));
plot([xx xx],[yy-2 yy+2],'-k','linewidth',2)
end
legend('Bouguer Anomaly', 'Faults', 'Location', 'best');
4 Comments
Mathieu NOE
on 27 Mar 2024
hello again
ok, so what is the problem with another data set ? the code would need to be tweaked or ?
What we are looking at are segments with a slope above a treshold value , if I understand correctly
do you fear that this approach would not be robust enough with another data set ?
More Answers (0)
See Also
Categories
Find more on Assumptions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!