In daily work, evaluating seismic activity often involves simple numerical analysis of the earthquake catalog regarding time, space, and intensity. Previously, I used Mapsis, which is very powerful, and Professor Li is great. However, the version of Mapsis I have is quite old, and it seems there is no new version available. I am not sure about the specifics, but the generated figures do not meet the expert review opinions, and I have to modify them in MapInfo every time, which is very painful. I don’t want to talk about it here. I am not very familiar with other functions, so I simply wrote some code in MATLAB with the attitude of just getting it done. The figures are ugly, but it saves time and avoids the hassle of constant modifications. Also, the safety evaluation report is very formatted; why isn’t there a formatting software tool? Maybe I just haven’t found one. If there is one, please recommend it, as spending money is always better than wasting people’s time.
The EQT format is as follows:
20010110225727 38.83 104.122.30 8 0
20010115191913 39.08 104.472.90 8 0
20010117030849 39.75 106.472.60 8 0
MATLAB displays as follows:
=== Earthquake Data Reading Process ===
Original data line 1: 20010110225727 38.83 104.122.30 8 0
Original data line 2: 20010115191913 39.08 104.472.90 8 0
Original data line 3: 20010117030849 39.75 106.472.60 8 0
Parsed result: Latitude 38.83, Longitude 104.12, Magnitude 2.3, Depth 8km
—————————-
Parsed result: Latitude 39.08, Longitude 104.47, Magnitude 2.9, Depth 8km
—————————-
Parsed result: Latitude 39.75, Longitude 106.47, Magnitude 2.6, Depth 8km
—————————-
=== Basic Statistical Information ===
Successfully read 1162 earthquake events
Average depth: 12.45km (1.00-45.00km)
Average magnitude: 2.50 (2.00-4.60)
Latitude range: 36.90°-40.00°
Longitude range: 104.12°-107.43°
Time range: January 10, 2001 – July 5, 2025
=== Distribution Statistics ===
Depth distribution:
1-4 km: 39 times (3.36%)
5-9 km: 562 times (48.36%)
10-14 km: 187 times (16.09%)
15-19 km: 144 times (12.39%)
20-24 km: 101 times (8.69%)
25-29 km: 76 times (6.54%)
30-34 km: 45 times (3.87%)
≥35 km: 8 times (0.69%)
Magnitude distribution:
M2.0-2.9: 965 times (83.05%)
M3.0-3.9: 169 times (14.54%)
M4.0-4.6: 28 times (2.41%)
=== Analysis Complete ===
Generated charts: 1. Depth distribution 2. 3D distribution 3. Profile analysis 4. Time profile
Statistical results saved: earthquake_statistics.txt
Charts saved as PNG files
MATLAB code is as follows:
clc; clear; close all;
%% 1. Read file
filename = '204611111.EQT';
fprintf('=== Earthquake Data Reading Process ===\n');
% Read all lines at once
lines = readlines(filename, 'EmptyLineRule', 'skip');
lines(cellfun(@isempty, lines)) = []; % Remove empty lines
nLines = length(lines);
% Preallocate arrays
maxEvents = nLines;
depths = zeros(maxEvents, 1); lats = zeros(maxEvents, 1);
lons = zeros(maxEvents, 1); mags = zeros(maxEvents, 1);
yrs = zeros(maxEvents, 1); months = zeros(maxEvents, 1);
days = zeros(maxEvents, 1); hours = zeros(maxEvents, 1);
minutes = zeros(maxEvents, 1); seconds = zeros(maxEvents, 1);
validCount = 0;
% Parse data
for i = 1:min(3, nLines) % Only show debugging information for the first 3 lines
fprintf('Original data line %d: %s\n', i, lines{i});
end
for i = 1:nLines
line = lines{i};
if length(line) < 35, continue; end % Skip incomplete lines
% Extract string fields
datetime_str = strtrim(line(1:15));
lat_str = strtrim(line(17:21));
lon_str = strtrim(line(23:28));
mag_str = strtrim(line(29:31));
depth_str = strtrim(line(34:35));
% Seconds 60 correction
if length(datetime_str) >= 14 && strcmp(datetime_str(13:14), '60')
datetime_str(13:14) = '59';
end
% Numeric conversion
[lat, lon, mag, depth] = deal(str2double(lat_str), str2double(lon_str), ...
str2double(mag_str), str2double(depth_str));
if all(~isnan([depth, lat, lon, mag])) && length(datetime_str) >= 14
validCount = validCount + 1;
% Store data
depths(validCount) = depth;
lats(validCount) = lat;
lons(validCount) = lon;
mags(validCount) = mag;
% Parse time
yrs(validCount) = str2double(datetime_str(1:4));
months(validCount) = str2double(datetime_str(5:6));
days(validCount) = str2double(datetime_str(7:8));
hours(validCount) = str2double(datetime_str(9:10));
minutes(validCount) = str2double(datetime_str(11:12));
seconds(validCount) = str2double(datetime_str(13:14));
% Display parsing results for the first 3 lines
if i <= 3
fprintf('Parsed result: Latitude %.2f, Longitude %.2f, Magnitude %.1f, Depth %.0fkm\n', ...
lat, lon, mag, depth);
fprintf('----------------------------\n');
end
end
end
% Truncate arrays to actual valid data count
depths = depths(1:validCount); lats = lats(1:validCount);
lons = lons(1:validCount); mags = mags(1:validCount);
yrs = yrs(1:validCount); months = months(1:validCount);
days = days(1:validCount); hours = hours(1:validCount);
minutes = minutes(1:validCount); seconds = seconds(1:validCount);
%% 2. Basic statistics
fprintf('\n=== Basic Statistical Information ===\n');
fprintf('Successfully read %d earthquake events\n', validCount);
fprintf('Average depth: %.2fkm (%.2f-%.2fkm)\n', mean(depths), min(depths), max(depths));
fprintf('Average magnitude: %.2f (%.2f-%.2f)\n', mean(mags), min(mags), max(mags));
fprintf('Latitude range: %.2f°-%.2f°\n', min(lats), max(lats));
fprintf('Longitude range: %.2f°-%.2f°\n', min(lons), max(lons));
if validCount > 0
minIdx = find(yrs == min(yrs), 1);
maxIdx = find(yrs == max(yrs), 1, 'last');
fprintf('Time range: %d year %d month %d day - %d year %d month %d day\n', ...
yrs(minIdx), months(minIdx), days(minIdx), ...
yrs(maxIdx), months(maxIdx), days(maxIdx));
end
%% 3. Depth and magnitude distribution statistics
% Depth distribution
depth_bins = [1,5,10,15,20,25,30,35,Inf];
depth_labels = {'1-4','5-9','10-14','15-19','20-24','25-29','30-34','≥35'};
depth_counts = histcounts(depths, depth_bins);
depth_pct = depth_counts / validCount * 100;
% Magnitude distribution
mag_bins = [2, 3, 4, 5];
mag_labels = {'M2.0-2.9', 'M3.0-3.9', 'M4.0-4.6'};
mag_counts = histcounts(mags, mag_bins);
mag_pct = mag_counts / validCount * 100;
% Unified display of statistical results
fprintf('\n=== Distribution Statistics ===\n');
fprintf('Depth distribution:\n');
for i = 1:length(depth_labels)
fprintf(' %s km: %d times (%.2f%%)\n', depth_labels{i}, depth_counts(i), depth_pct(i));
end
fprintf('Magnitude distribution:\n');
for i = 1:length(mag_labels)
fprintf(' %s: %d times (%.2f%%)\n', mag_labels{i}, mag_counts(i), mag_pct(i));
end
%% 4. Chart parameters
colors = [0.2 0.8 0.2; 0.9 0.9 0.2; 0.8 0.2 0.2]; % Green, Yellow, Red
sizes = [30, 60, 90];
magBin = discretize(mags, mag_bins);
%% 5. Chart 1: Depth distribution percentage bar chart
figure('Position', [100 100 900 600], 'Color', 'w');
bar(depth_pct, 'FaceColor', [0.27 0.62 0.81], 'EdgeColor', [0.1 0.3 0.5], ...
'LineWidth', 1.2, 'FaceAlpha', 0.8);
set(gca, 'FontSize', 12, 'FontWeight', 'bold', 'XTickLabel', depth_labels, ...
'GridAlpha', 0.3, 'LineWidth', 1);
grid on
xlabel('Source Depth (km)', 'FontSize', 14, 'FontWeight', 'bold');
ylabel('Percentage (%)', 'FontSize', 14, 'FontWeight', 'bold');
title(sprintf('Source Depth Distribution Statistics Total Earthquake Count: %d', validCount), ...
'FontSize', 16, 'FontWeight', 'bold');
% Add value annotations
for i = 1:length(depth_counts)
text(i, depth_pct(i) + max(depth_pct)*0.07, sprintf('%d times', depth_counts(i)), ...
'Horizontal', 'center', 'FontSize', 11, 'FontWeight', 'bold');
text(i, depth_pct(i) + max(depth_pct)*0.03, sprintf('%.2f%%', depth_pct(i)), ...
'Horizontal', 'center', 'FontSize', 11, 'FontWeight', 'bold');
end
ylim([0 max(depth_pct)*1.15]);
%% 6. Chart 2: 3D source distribution chart
figure('Position', [200 200 1000 700], 'Color', 'w');
hold on
scatterHandles = gobjects(1, 3);
for k = 1:3
idx = (magBin == k);
if any(idx)
scatterHandles(k) = scatter3(lons(idx), lats(idx), -depths(idx), ...
sizes(k), colors(k,:), 'filled', ...
'MarkerEdgeColor', 'k', 'LineWidth', 0.5);
end
end
xlabel('Longitude (°)'); ylabel('Latitude (°)'); zlabel('Depth (km)');
title('3D Source Distribution Chart'); grid on; box on; view(-35, 25);
legend(scatterHandles, mag_labels, 'Location', 'best');
%% 7. Chart 3: Profile analysis chart
figure('Position', [300 300 1200 800], 'Color', 'w');
% Calculate diagonal profile
start_point = [min(lons), max(lats)]; % Northwest
end_point = [max(lons), min(lats)]; % Southeast
[diag_length_km, ~] = distance(start_point(2), start_point(1), end_point(2), end_point(1));
diag_length_km = deg2km(diag_length_km);
% Calculate projected distances
diag_vec = end_point - start_point;
diag_unit = diag_vec / norm(diag_vec);
proj_dists = zeros(validCount, 1);
for i = 1:validCount
point_vec = [lons(i), lats(i)] - start_point;
proj_len = dot(point_vec, diag_unit);
proj_dists(i) = max(0, min(norm(diag_vec), proj_len));
end
proj_dists_km = (proj_dists / norm(diag_vec)) * diag_length_km;
% Upper subplot: Depth profile
subplot(2,1,1); hold on
for k = 1:3
idx = (magBin == k);
if any(idx)
scatter(proj_dists_km(idx), -depths(idx), sizes(k), colors(k,:), ...
'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 0.5);
end
end
xlabel('Distance along profile line (km)'); ylabel('Depth (km)');
title('Source Depth Distribution along Profile Line'); grid on; box on;
xlim([0, diag_length_km]);
% Add profile line schematic
hold on
plot([0, diag_length_km], [min(-depths)-5, min(-depths)-5], 'r-', 'LineWidth', 3, ...
'Color', 'red');
plot(0, min(-depths)-5, '^', 'MarkerSize', 8, ...
'MarkerFaceColor', 'blue', 'MarkerEdgeColor', 'k');
plot(diag_length_km, min(-depths)-5, 'v', 'MarkerSize', 8, ...
'MarkerFaceColor', 'red', 'MarkerEdgeColor', 'k');
% Lower subplot: Plane distribution
subplot(2,1,2); hold on
plotHandles = gobjects(1, 3);
for k = 1:3
idx = (magBin == k);
if any(idx)
plotHandles(k) = scatter(lons(idx), lats(idx), sizes(k), colors(k,:), ...
'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 0.5);
end
end
% Draw profile line
diag_line = plot([start_point(1), end_point(1)], [start_point(2), end_point(2)], ...
'r-', 'LineWidth', 3);
start_marker = plot(start_point(1), start_point(2), '^', 'MarkerSize', 10, ...
'MarkerFaceColor', 'blue', 'MarkerEdgeColor', 'k');
end_marker = plot(end_point(1), end_point(2), 'v', 'MarkerSize', 10, ...
'MarkerFaceColor', 'red', 'MarkerEdgeColor', 'k');
xlabel('Longitude (°)'); ylabel('Latitude (°)');
title('Source Plane Distribution Chart and Profile Line'); grid on; box on;
legend([plotHandles, diag_line, start_marker, end_marker], ...
[mag_labels, {'Profile Line', 'Start Point', 'End Point'}], 'Location', 'best');
%% 8. Chart 4: Time-depth profile chart (if time data is available)
if validCount > 0 && any(yrs > 0)
figure('Position', [400 400 1200 500], 'Color', 'w');
hold on
timeHandles = gobjects(1, 3);
for k = 1:3
idx = (magBin == k);
if any(idx)
timeHandles(k) = scatter(yrs(idx), -depths(idx), sizes(k), colors(k,:), ...
'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 0.5);
end
end
xlabel('Year'); ylabel('Depth (km)'); title('Source Time-Depth Profile Chart');
grid on; box on; legend(timeHandles, mag_labels, 'Location', 'best');
unique_yrs = unique(yrs);
if length(unique_yrs) > 10
xticks(unique_yrs(1:ceil(length(unique_yrs)/10):end));
else
xticks(unique_yrs);
end
xtickangle(45);
end
%% 9. Save results
fprintf('\n=== Analysis Complete ===\n');
fprintf('Generated charts: 1. Depth distribution 2. 3D distribution 3. Profile analysis');
if validCount > 0 && any(yrs > 0), fprintf(' 4. Time profile'); end
fprintf('\n');
% Save statistical results
stats_file = 'earthquake_statistics.txt';
fid = fopen(stats_file, 'w');
fprintf(fid, 'Earthquake Statistical Analysis Report\nGenerated Time: %s\nData File: %s\n\n', datestr(now), filename);
fprintf(fid, 'Basic Statistics:\nTotal Events: %d\n', validCount);
fprintf(fid, 'Depth Range: %.1f-%.1fkm (Average: %.2fkm)\n', min(depths), max(depths), mean(depths));
fprintf(fid, 'Magnitude Range: %.1f-%.1f (Average: %.2f)\n\n', min(mags), max(mags), mean(mags));
fprintf(fid, 'Depth Distribution:\n');
for i = 1:length(depth_labels)
fprintf(fid, '%s km: %d times (%.2f%%)\n', depth_labels{i}, depth_counts(i), depth_pct(i));
end
fprintf(fid, '\nMagnitude Distribution:\n');
for i = 1:length(mag_labels)
fprintf(fid, '%s: %d times (%.2f%%)\n', mag_labels{i}, mag_counts(i), mag_pct(i));
end
fclose(fid);
fprintf('Statistical results saved: %s\n', stats_file);
% Save charts
try {
saveas(1, 'depth_distribution.png'); saveas(2, '3d_distribution.png');
saveas(3, 'spatial_analysis.png');
if validCount > 0 && any(yrs > 0), saveas(4, 'time_depth_profile.png'); end
fprintf('Charts saved as PNG files\n');
} catch {
fprintf('Chart saving failed\n');
}