plot earthquake locations data on a geotiff file

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)])

Answer 1

[Dyuman Joshi]

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.

[Anitha Limann]

Hello,

I did try these functions but could not get it to plot. Attached here is the data folder. Thank you.

[Dyuman Joshi]

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)])

[Anitha Limann]

Thank you a lot. This works.

I do have one question though. Is it possible to color code the earthquake location based on their depths?

Data = load('data/A_1.txt','-ascii'); % earthquake locations
lon = Data(:,1); 
lon = 360+lon; % longitude
lat = Data(:,2); % latitude
depth = Data(:,3); % Depth

When I use "plot" command I created a colormap to color code each earthquake based on its depth. But then I could not add a basemap to that plot.

Now I cannot use my colormap with this "geoshow" command. Sorry, I am not quite familiar with how "geoshow" works.

[Dyuman Joshi]

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'

[Anitha Limann]

Thank you for undersatnding and helping me with this.

Code below is the colormap I created to apply here. I tried to change it as you explained earlier for geoshow. But I am not sure where I m doing it wrong.

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);

% % 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), 'Color', c(ii,:),'LineStyle', 'none', 'Marker', '.', 'MarkerSize', 5);
    
end

[Dyuman Joshi]

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]

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

[Dyuman Joshi]

You are welcome!