plot earthquake locations data on a geotiff file

Question

Anitha Limann 2023-10-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2034519-plot-earthquake-locations-data-on-a-geotiff-file

评论： Dyuman Joshi 2023-10-20

Hello,

I want to plot locations of earthquakes on a base map. I do have a geotiff file that i can use to create a base map but using following code i was unable to map the points on the base map. Please let me know the error in this code.

Data = load('Area_1.txt','-ascii'); % earthquake locations
lon = Data(:,1); lon = 360+lon; % longitude
lat = Data(:,2); % latitude
nv1 = importdata("volcanic.xyz"); % spreading ridge
lon1 = nv1(:,1); % longitude
lat1 = nv1(:,2); % latitude
figure(1); clf; axis equal; hold on;
    
    basemap = readgeoraster("exportImage.tiff");
    imshow(basemap); 
plot(lon1, lat1,'.'); % spreading ridge
scatter(lon, lat, 15, c,'filled','o','LineWidth', 1, 'MarkerFaceAlpha', 0.9); % earthquake locations
ylabel('latitude (deg)');   xlabel('longitude (deg)');
title(['CMT Depth distribution ', Name])
xlim([min(lon) max(lon)]);    ylim([min(lat) max(lat)])

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Dyuman Joshi 2023-10-17

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2034519-plot-earthquake-locations-data-on-a-geotiff-file#answer_1334594

You should be using geoplot, geoscatter and mapshow/geoshow instead of plot(), scatter() and imshow().

If, even after making the above mentioned changes, the code is not working or the output is not as desired, please attach the data files you are working with.

8 个评论
显示 6更早的评论隐藏 6更早的评论

Dyuman Joshi 2023-10-17

编辑：Dyuman Joshi 2023-10-17

在 MATLAB Online 中打开

data.zip

I failed to recall that geoplot/geoscatter do not work with geoshow while suggesting them above. I have plotted all the data using geoshow().

@Anitha Limann, please check out the output below and let me know if it works or not.

%Unzip and get the file names

y = unzip('data.zip')

y = 1×4 cell array

{'data/'} {'data/A_1.txt'} {'data/volcanic.xyz'} {'data/exportImage.tiff'}

Data = load('data/A_1.txt','-ascii'); % earthquake locations

lon = Data(:,1);

lon = 360+lon; % longitude

lat = Data(:,2); % latitude

nv1 = importdata("data/volcanic.xyz"); % spreading ridge

lon1 = nv1(:,1); % longitude

lat1 = nv1(:,2); % latitude

%Add the image as base

geoshow("data/exportImage.tiff");

hold on

%Add the data

geoshow(lat1, lon1, 'DisplayType', 'line', 'Color', [0 1 0], 'LineWidth', 1.5) % spreading ridge

geoshow(lat, lon, 'DisplayType', 'line', 'Color', [0 1 1], 'LineStyle', 'none', 'Marker', '.', 'MarkerSize', 10)

ylabel('latitude (deg)');

xlabel('longitude (deg)');

title(['CMT Depth distribution'])

%Sample legend

%modify as necessary

legend({'','Ridge','Volcanos'}, 'location', 'southeast')

%geolimits([min(lat) max(lat)], [min(lon) max(lon)])

Dyuman Joshi 2023-10-17

在 MATLAB Online 中打开

data.zip

Don't be sorry, it's not easy to work with as well.

It seems a functionality of geoshow() is not working properly (which I have specified at the bottom), so I have used a for loop as a workaround. However, going by this method the legends will be according to the depth.

%Unzip and get the file names

y = unzip('data.zip');

Data = load('data/A_1.txt','-ascii'); % earthquake locations

lon = Data(:,1);

lon = 360+lon; % longitude

lat = Data(:,2); % latitude

depth = Data(:,3);

%Find the unique depth values

d = unique(depth);

n = numel(d);

nv1 = importdata("data/volcanic.xyz"); % spreading ridge

lon1 = nv1(:,1); % longitude

lat1 = nv1(:,2); % latitude

%Add the image as base

geoshow("data/exportImage.tiff");

hold on

%colormap

%random map chosen for example

map = spring(n);

%Add the data

geoshow(lat1, lon1, 'DisplayType', 'line', 'Color', [0 1 1], ...

'LineWidth', 1.5) % spreading ridge

for k=1:n

idx=depth==d(k);

geoshow(lat(idx), lon(idx), 'DisplayType', 'line', 'Color', map(k,:), ...

'LineStyle', 'none', 'Marker', '.', 'MarkerSize', 10)

end

ylabel('latitude (deg)');

xlabel('longitude (deg)');

title(['CMT Depth distribution'])

%Define string for legend

%Add unit to the depth accordingly

str = [""; "Ridge"; compose("Volcanoes at depth = %d units", d)]

str = 5×1 string array

"" "Ridge" "Volcanoes at depth = 12 units" "Volcanoes at depth = 15 units" "Volcanoes at depth = 22 units"

legend(str, 'location', 'southeast')

Weird behaviour - Initializing geoshow with multipoint (similar to scatter) does not seem to be working. It still registers the data as a line (similar to plot)

Reference - https://in.mathworks.com/help/map/ref/geoshow.html?s_tid=doc_ta#d126e68241

figure

geoshow("data/exportImage.tiff");

s = geoshow(lat, lon, 'DisplayType', 'multipoint');

s.Type

ans = 'line'

Dyuman Joshi 2023-10-18

在 MATLAB Online 中打开

data.zip

Here you go -

%Unzip and get the file names

y = unzip('data.zip');

Data = load('data/A_1.txt','-ascii'); % earthquake locations

lon = Data(:,1);

lon = 360+lon; % longitude

lat = Data(:,2); % latitude

depth = Data(:,3);

%Added snippet of code

depth_levels = 10:5:30;

c1 = [244 163 234]./256; c2 = [180 60 161]./256;

c3 = [106 17 179]./256; c4 = [0 0 0]./256;

colors = [c1;c2;c3;c4];

depth_levels_centers = (depth_levels(1:end-1)+depth_levels(2:end))/2;

c = zeros(size(Data,1),1);

nv1 = importdata("data/volcanic.xyz"); % spreading ridge

lon1 = nv1(:,1); % longitude

lat1 = nv1(:,2); % latitude

%% Add the image as base

%Turning the handlevisibility off as to not show the legend for the image

geoshow("data/exportImage.tiff", 'HandleVisibility', 'off');

%Display the data

geoshow(lat1, lon1, 'DisplayType', 'line', 'Color', [0 1 1], ...

'LineWidth', 1.5, 'DisplayName', 'Ridge') % spreading ridge

hold on

%% for each depth range ...

for ii = 1:numel(depth_levels)-1

% get a logical index of the data points in that range

idx = depth > depth_levels(ii) & depth <= depth_levels(ii+1);

c(idx,:) = depth_levels_centers(ii);

geoshow(lat(idx), lon(idx), 'DisplayType', 'line', 'Color', colors(ii,:), ...

'LineStyle', 'none', 'Marker', '.', 'MarkerSize', 10, ...

'DisplayName', sprintf('Volcanoes in the depth range = %d - %d units', depth_levels(ii:ii+1)));

end

%Text and labels

ylabel('latitude (deg)');

xlabel('longitude (deg)');

title(['CMT Depth distribution'])

legend('Location', 'southeast')

Anitha Limann 2023-10-20

Thank you so much for your kind support. You have been extremely helpful to me. Thanks again.

Dyuman Joshi 2023-10-20

You are welcome!

请先登录，再进行评论。

plot earthquake locations data on a geotiff file

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

8 个评论
显示 6更早的评论隐藏 6更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

plot earthquake locations data on a geotiff file

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

8 个评论 显示 6更早的评论隐藏 6更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

WeChat

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论