I for got the files
Need to put symbols on curve
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Moustafa Abedel Fattah
le 26 Mar 2024
Commenté : Mathieu NOE
le 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');
1 commentaire
Réponse acceptée
Mathieu NOE
le 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 commentaires
Mathieu NOE
le 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 ?
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Assumptions dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!