フィルターのクリア

Need to put symbols on curve

10 ビュー (過去 30 日間)
Moustafa Abedel Fattah
Moustafa Abedel Fattah 2024 年 3 月 26 日
コメント済み: Mathieu NOE 2024 年 3 月 27 日
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:');
Detected Faults:
disp('Location Angle (degrees) Type');
Location Angle (degrees) Type
for i = 1:numel(fault_locations)
disp([num2str(xc(fault_locations(i))), ' ', num2str(fault_angles(i)), ' ', fault_types{i}]);
end
32.5335 87.5575 Normal 34.1602 87.4417 Normal 35.7869 87.3003 Normal 37.4135 86.9421 Normal 39.0402 86.7289 Normal 40.6669 86.5637 Normal 42.2936 85.9853 Normal 43.9202 85.5759 Normal 45.5469 85.5344 Normal 47.1736 85.3 Normal 48.8003 85.1095 Normal 50.4269 85.0995 Normal 52.0536 84.9447 Normal 53.6803 84.8022 Normal 55.307 84.8042 Normal 56.9336 84.6707 Normal 58.5603 84.5332 Normal 60.187 84.5309 Normal 61.8137 84.0933 Normal 63.4403 83.5875 Normal 65.067 83.5192 Normal 66.6937 83.2219 Normal 68.3204 82.8312 Normal 69.947 82.7367 Normal 71.5737 82.5562 Normal 73.2004 82.2813 Normal 74.8271 82.1869 Normal 76.4537 81.9931 Normal 78.0804 81.6392 Normal 79.7071 81.4811 Normal 81.3338 81.6211 Normal 82.9604 81.935 Normal 84.5871 82.1116 Normal 86.2138 82.3994 Normal 87.8405 83.2564 Normal 89.4671 83.8262 Normal 91.0938 83.8981 Normal 92.7205 84.2241 Normal 94.3472 84.4756 Normal 95.9738 84.4911 Normal 97.6005 84.634 Normal 99.2272 84.7591 Normal 100.8539 84.7591 Normal 102.4805 84.8852 Normal 104.1072 85.0095 Normal 105.7339 85.0079 Normal 107.3606 85.0745 Normal 108.9872 85.1483 Normal 110.6139 85.1542 Normal 112.2406 85.1973 Normal 113.8673 85.2511 Normal 115.4939 85.2608 Normal 117.1206 85.2625 Normal 118.7473 85.2651 Normal 120.374 85.2644 Normal 122.0006 85.2911 Normal 123.6273 85.3379 Normal 125.254 85.3566 Normal 126.8807 85.4039 Normal 128.5073 85.5032 Normal 130.134 85.5545 Normal 131.7607 85.5996 Normal 133.3874 85.72 Normal 135.014 85.7957 Normal 136.6407 85.8559 Normal 138.2674 86.0763 Normal 139.8941 86.2382 Normal 141.5207 86.4028 Normal 143.1474 87.3928 Normal 250.508 87.7356 Normal 252.1346 87.722 Normal 253.7613 87.5996 Normal 255.388 87.4512 Normal 257.0147 87.4251 Normal 258.6413 87.3466 Normal 260.268 87.2373 Normal 261.8947 87.2064 Normal 263.5214 87.1593 Normal 265.148 87.0823 Normal 266.7747 87.0521 Normal 268.4014 87.0475 Normal 270.0281 87.0386 Normal 271.6547 87.0335 Normal 273.2814 87.0614 Normal 274.9081 87.1328 Normal 276.5348 87.1752 Normal 278.1615 87.2286 Normal 279.7881 87.4138 Normal 281.4148 87.5449 Normal 283.0415 87.6241 Normal 310.695 86.8596 Normal 312.3216 86.5014 Normal 313.9483 86.029 Normal 315.575 85.1678 Normal 317.2017 84.7776 Normal 318.8283 84.9774 Normal 320.455 85.4346 Normal 322.0817 85.6909 Normal 323.7084 85.9865 Normal 325.335 0
% 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 件のコメント
Moustafa Abedel Fattah
Moustafa Abedel Fattah 2024 年 3 月 26 日
移動済み: Voss 2024 年 3 月 26 日

サインインしてコメントする。

採用された回答

Mathieu NOE
Mathieu NOE 2024 年 3 月 26 日
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 件のコメント
Moustafa Abedel Fattah
Moustafa Abedel Fattah 2024 年 3 月 26 日
Thank yo I mean fha we needint the red symboles its enough the five-vertical line as indicated for the location of fault. but I need to located it automatical if I load another data of (xc, gBt)
Thanks for your accepted answer
Mathieu NOE
Mathieu NOE 2024 年 3 月 27 日
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 ?

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeAssumptions についてさらに検索

製品


リリース

R2014b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by