How plot line in geographic map with initial coordinates and azimuth?

Hello!

I would like to know how can I plot a line extending for 3000 km over a geographic map on MATLAB. I have the latitude and longitude of the initial position, together of the back-azimuth to indicate the direction. Then, I need to plot a line starting in that position. I know some commands (e.g.geoplot) which they plot line between two geographical coordinates, but I dont know using bach-azimuth angle.

I will be grateful for any idea.

best,

Guilherme de Melo

Answer 1

[Cris LaPierre]

I think geoplot is the correct function. You just need to calculate the end coordinates. If you have the initial location, arclength, and azimuth, you can use the reckon function (included in the Mapping Toolbox).

[Guilherme de Melo]

Hi Cris,

It worked using the reckon, and I got to create a variable with all the final latitude and longitudes coordinates. Now, I am trying to plot all the lines in map, however, it is not working. It is plotting a single line.

figure(12)
k=0;
for i=1:N
    k=k+1;
    geoplot([lat_mean_st123_array_plot(k) lat_final_line_point(k)],[long_mean_st123_array_plot(k) long_final_line_point(k)],'-')
end
geobasemap satellite

lat_mean_st123_array is the latitude where is located the array, and long_mean_st123_array is the longitude. Each one have a single coordinate. lat_final_line_point and long_final_line_point have 17 coordinates each one, which refer the final position of the 17 lines.

Do you know what I need to change? I will be grateful.

thanks,

Guilherme

[Cris LaPierre]

You need to include a hold on to plot multiple lines on the same axes. Be sure to include a hold off outside your loop. i've also simplified the code a little.

figure(12)
for k=1:N
    geoplot([lat_mean_st123_array_plot(k) lat_final_line_point(k)],[long_mean_st123_array_plot(k) long_final_line_point(k)],'-')
    hold on
end
hold off
geobasemap satellite

[Guilherme de Melo]

Hi Cris,

Thank you very much. It worked very well. I got to do also using the "plot" function:

figure(8)
%text(lat_new_epicenter,lon_new_epicenter,labels,'Color','red')
gx = geoaxes;
plot(gx,lat,long,lat_globalnet_epicenter_plotline,long_globalnet_epicenter_plotline,"k*",'Color','white')
hold on
plot(gx,lat_mean_line,long_mean_line,'Color','red','LineWidth',2)
geobasemap satellite

Do you know how can I apply different color to each line? For example, applying some pallete color in total lines? The number of the lines varies depending of the parameters fixed at the start of my signal analysis. So, the total number of lines can be sometimes small and large.

[Cris LaPierre]

Use the colororder function to set a predetermined list of colors your plots cycle through.

You can also extract a specific number of unique colors from a particular colormap. For example, parula

C = parula(5)

C = 5×3

0.2422 0.1504 0.6603 0.2021 0.4788 0.9911 0.0704 0.7457 0.7258 0.7858 0.7573 0.1598 0.9769 0.9839 0.0805

[Guilherme de Melo]

Hi Cris,

I tried using the following:

C = parula(i);
plot(gx,lat,long,lat_globalnet_epicenter_plotline,long_globalnet_epicenter_plotline,"k*",'Color',[C])

"i" have the total number of segments (19), which is the same total number of lat, log, lat_globalnet_epicenter_plotline, long_globalnet_epicenter_plotline, then the total 19 lines should be plotted in different color each one. However, it is showing the following error:

Error using plot

Invalid RGB triplet. Specify a three-element vector of values between 0 and 1.

Error in source_location_updated (line 875)

plot(gx,lat,long,lat_globalnet_epicenter_plotline,long_globalnet_epicenter_plotline,"k*",'Color',[C])

[Cris LaPierre]

Just a note that we do not see what your plot looks like, but you must either specify a color for each series separately in the plot command, or use the colororder function to specify a color order for the axes.

I would suggest the following. This will plot each data series using the colors in C. Each column of data is a data series.

C = parula(i);
colororder(C);
plot(gx,lat,long,lat_globalnet_epicenter_plotline,long_globalnet_epicenter_plotline,"k*") 

[Guilherme de Melo]

Thank you. It worked well. Image is attached. Here was the code:

figure(8)
%text(lat_new_epicenter,lon_new_epicenter,labels,'Color','red')
gx = geoaxes;
C = parula(N);
colororder(C);
plot(gx,lat,long)
hold on
plot(gx,lat_globalnet_epicenter_plotline,long_globalnet_epicenter_plotline,"k*",'Color','white')
hold on
plot(gx,lat_mean_line,long_mean_line,'Color','red','LineWidth',2)
colorbar
geobasemap satellite
geolimits([str2num(coord_lat_range_tracklines_map)],[str2num(coord_long_range_tracklines_map)])
print(gcf,'Fig8_baz_lines_map.png','-dpng','-r300');

Do you works with xcorr command of Matlab?

[Cris LaPierre]

I have used it, but I do not work on that team.

[Guilherme de Melo]

Hi Cris,

I got to plot using the reckon and it worked well. He is the code:

distance_line=track_lines_distance; %km
distUnits = 'km';
arclen_array_line = rad2deg((distance_line)/earthRadius(distUnits));
[lat_final_line_points,long_final_line_points] = reckon(lat_mean_st123_array, long_mean_st123_array,arclen_array_line,baz_lines)

All lines have the same extension provided on "track_lines_distance" variable. For example, I defined 1000 km and they were plotted well. However, do you know how I can determine to it decrease each line for a specific value? I mean, "back_line" has 20 azimuths, generating 20 lines of 1000 km extension. But, I would like that it was 1000, 1990, 1980, 1970....suppossing that I create a variable "remove_line_extension=10;"

Do you know how can I do that? I will be very grateful.

best,

Guilherme

[Cris LaPierre]

What you are describing doesn't sound like standard practice, which means it may be possible, but likely requires a custom solution. There likely is not a helper function or input argument for that specific need.

It would be easier to try out if you could share a working example of what you have. You can add executable code to your post using the Code option in the ribbon, and can attach files using the paperclip icon.

[Guilherme de Melo]

Hi,

Ok. You can run the example code attached, and you will see a map and three lines. You can observe the lines with the same final position, approximately.

clc; clear; close all;

track_lines_distance=1000; %line extension in km
lat_mean_st123_array=-15; %latitude
long_mean_st123_array=-10; %longitude
baz_lines=[268;270;271]; %example of back-azimuths

distance_line=track_lines_distance; %km
distUnits = 'km';
arclen_array_line = rad2deg((distance_line)/earthRadius(distUnits));
[lat_final_line_point,long_final_line_point] = reckon(lat_mean_st123_array, long_mean_st123_array,arclen_array_line,baz_lines);

lat_mean_st123_array_plot=[-15,-15,-15];
long_mean_st123_array_plot=[-10,-10,-10];


lat=[lat_mean_st123_array_plot; transpose(lat_final_line_point)];
long=[long_mean_st123_array_plot; transpose(long_final_line_point)];

figure(1)
gx = geoaxes;
plot(gx,lat,long)
geobasemap grayterrain

[Guilherme de Melo]

I mean, horizontal position (longitude).

[Cris LaPierre]

Simplest way I can think of is to call reckon multiple times with the desired distance.

track_lines_distance=1000; %line extension in km
lat_mean_st123_array=-15; %latitude
long_mean_st123_array=-10; %longitude
baz_lines=[268;270;271]; %example of back-azimuths
remove_line_extension=10;

distUnits = 'km';

lat_final_line_point=[];
long_final_line_point=[];
for b = 1:length(baz_lines)
    distance_line=track_lines_distance - (b-1)*remove_line_extension; %km
    arclen_array_line = rad2deg((distance_line)/earthRadius(distUnits));
    [lat_final_line_point(1,b),long_final_line_point(1,b)] = reckon(lat_mean_st123_array, long_mean_st123_array,arclen_array_line,baz_lines(b));
end

lat_mean_st123_array_plot=[-15,-15,-15];
long_mean_st123_array_plot=[-10,-10,-10];

lat=[lat_mean_st123_array_plot; lat_final_line_point];
long=[long_mean_st123_array_plot; long_final_line_point];

figure(1)
gx = geoaxes;
plot(gx,lat,long)

[Guilherme de Melo]

Hi

Thank you very much. The lines decreased the extension parallel to the azimuths change in the array. :-))

There is another array with time referent to each line. For example:

First line - 1 s

Second line - 2 s

Third line - 3 s

....

I tried to color the lines using a palette color with the lines changing of color from the azimuth and line extension variation. It means the first line starts in dark blue and last with dark red. If you take a look at the attached picture, it is not in agreement. Do you know why? I show you in follow the plotted figure code:

figure(8)

gx = geoaxes;
plot(gx,lat,long,'LineWidth',3.5)

colormap(jet); %preparing palete color
t = colorbar('Ticks',time_vector);
set(get(t,'ylabel'),'String','Time-window order, s','fontsize',14); %defining the label
caxis([min(time_vector) max(time_vector)]); %defining the color axis from the time windows number

geobasemap grayterrain

[Cris LaPierre]

colormap sets the color palette, but does not set the color order used for plotting. For that, use colororder. By default, there are 7 colors. Once you cycle through all the colors, it restarts at the 1st color. That means the 8th line is plotted in the same color as the 1st line. Increase the colors in colororder to ensure each line is plotted in a unique color.

[Guilherme de Melo]

Hi. I updated my R2020 to R2023 to get the colororder option, but it really do repetition after 8th color. is it possible to include a pallete color inside the coloroder? or how could I divide the pallete of colororder to the exact number of lines? Thanks

[DGM]

The default colororder is lines(), which only has 7 unique colors. Beyond that, they're just repeated.

get(groot,'DefaultAxesColorOrder')

ans = 7×3

0 0.4470 0.7410 0.8500 0.3250 0.0980 0.9290 0.6940 0.1250 0.4940 0.1840 0.5560 0.4660 0.6740 0.1880 0.3010 0.7450 0.9330 0.6350 0.0780 0.1840

If you want to have more unique colors, use a different color table generation tool (e.g. jet(), colorcube(), etc) or create a custom color Mx3 color table. There are also other tools on the File Exchange.

https://www.mathworks.com/matlabcentral/fileexchange/45208-colorbrewer-attractive-and-distinctive-colormaps

https://www.mathworks.com/matlabcentral/fileexchange/62729-matplotlib-perceptually-uniform-colormaps

https://www.mathworks.com/matlabcentral/fileexchange/29702-generate-maximally-perceptually-distinct-colors

https://www.mathworks.com/matlabcentral/fileexchange/70215-maximally-distinct-color-generator

If you want to use jet() (or other tools) to set the colororder, specify the number of colors that you actually need when you call the function.

% some fake data
numlines = 10;
x = 0:360;
y = sind(x).*linspace(0,1,numlines).';

% plot the lines
% set the colororder using the default CT length
% here, jet() returns 256 colors, but we're only using the first 10
plot(x,y,'linewidth',2)
colororder(jet)

figure
% plot the lines
% set the colororder using an appropriate length
% here, jet returns 10 colors, so we get the full sweep
plot(x,y,'linewidth',2)
colororder(jet(numlines))

[Guilherme de Melo]

Hi!

Thank you. The color lines were plotted good with base in the pallete jet. I used this:

total_elements_time_vector=numel(time_vector);
colororder(jet(total_elements_time_vector));

However, the pallete colorbar was not plotted in jet color, please check the atached figure. I tried using this command:

t = colorbar('Ticks',time_vector);
set(get(t,'ylabel'),'String','Time-window order, s','fontsize',14); %defining the label

Do you know what is the error?

[Cris LaPierre]

You shouldn't need to update. colororder was introduced in R2019b. However, there is no harm in doing so if you can.

The colorbar color comes from the colormap. that is not tied to colororder, so unless you explicitely change it, it will be using the default. You need to set the colormap to jet, too.

track_lines_distance=1000; %line extension in km
lat_mean_st123_array=-15; %latitude
long_mean_st123_array=-10; %longitude
baz_lines=268:280; %example of back-azimuths

distance_line=track_lines_distance; %km
distUnits = 'km';
arclen_array_line = rad2deg((distance_line)/earthRadius(distUnits));
[lat_final_line_point,long_final_line_point] = reckon(lat_mean_st123_array, long_mean_st123_array,arclen_array_line,baz_lines);

lat_mean_st123_array_plot=repmat(lat_mean_st123_array,size(baz_lines));
long_mean_st123_array_plot=repmat(long_mean_st123_array,size(baz_lines));


lat=[lat_mean_st123_array_plot; lat_final_line_point];
long=[long_mean_st123_array_plot; long_final_line_point];

% without setting colormap
figure(1)
colororder(jet(numel(baz_lines)))
gx = geoaxes;
plot(gx,lat,long)
geobasemap grayterrain
colorbar

% setting colormap
figure(2)
colororder(jet(numel(baz_lines)))
gx = geoaxes;
plot(gx,lat,long)
geobasemap grayterrain
colormap jet
colorbar

[Guilherme de Melo]

Hi Cris, thank you very much. I updated to R2023 because there are some color pallete of colororder available only on R2023 (https://www.mathworks.com/help/matlab/ref/colororder.html). However, my R2023 was stuck when I try to run the code. I kept the R2020 version.

[Guilherme de Melo]

Do you know how I can plot the same lines but using great circle lines projection? I mean, there are two coordinates where if you look a spherical global map, the direct line is always on sea. However, if we plot the two coordinates on 'plot' function of Matlab its is like a direct line crossing continents. Did you understand me?

[Cris LaPierre]

I believe that is called the Gnomonic Projection. See here: https://www.mathworks.com/help/map/summary-and-guide-to-projections.html#mw_39de9c72-1c2d-4011-ba0c-805090d20c77

[Guilherme de Melo]

Hi Cris, I got to plot map with gnomonic projection using the M_Map package for MatLab(https://www.eoas.ubc.ca/~rich/map.html).

I have another questions. Is there a way for us to use the reckon to project the lines until another specific array of coordinates? I mean, the horizontal thicker black line shown on the map (see the figure). I have all the coordinates of the line. Then, I would like to try to calculate the horizontal length (km) over the thicker black line between the further east and west projected lines. For example, taking a look at map scale the length should be ~70-74 km, but I am not sure about the exact length.

thanks,

[Cris LaPierre]

No, Reckon does not find the intercept of 2 lines. You will find this forum full of posts on how to do that.

One brute force approach if you have the coordinates of the black line might be to use the distance function to calculate what the distance to each point is, and then coding up a solution for picking the best point.

[Guilherme de Melo]

I took a look in several posts but seems that they are more curves or mathematical lines and not geographic coordinates. I am not sure if it will be possible on MATLAB. :-(

[Cris LaPierre]

If by not possible, you mean there is no built-in functionality for clippilng the trajectory when it hits the black line, you are correct. However, that does not mean it is not possible. It just means you will need to code an algorithm yourself that uses the built in capabilities of MATLAB to determine what distance to use.

My recommendation is to start a new question just about this. Be sure to include all the information you have. Share code that others can run that will exactly recreate your plot so that they can play around with it and test ideas. It is much harder to come up with a solution that will work if you only provide an image.

[Guilherme de Melo]

Hi Cris,

Do you know why the symbols are not colored from the colorbar which I ordered? back_azm_ev2_error has the azimuth errors. Then, I am trying to plot star symbols of the azimuth versus time, but with the symbols colored by the error. However, the symbols are showing by the single color like initially. thanks in advance

figure()
total_elements_baz_lines_ev2=numel(back_azm_ev2_error);
colororder(jet(total_elements_baz_lines_ev2)); %preparing palete color
plot(time_vector_ev2,baz_lines_ev2,'*')
colormap jet
t = colorbar('Ticks',[min(back_azm_ev2_error) max(back_azm_ev2_error)]);
set(get(t,'ylabel'),'String','Error, degree','fontsize',14);
caxis([0 10]);

[Cris LaPierre]

Most likely because your errors are all plotted as part of the same series. You'd have to share your variables if you want an answer specific to your data.

x1 = 1:5;
x2 = x1+1;
y = 1:5;

% all part of the same series
plot(x1,y,'*')
hold on
% Each point plotted its own series
for d = 1:length(y)
    plot(x2(d),y(d),'*')
end
hold off

As far as I can tell, there is no way in MATLAB to plot a vector of numbers as individual series (each with a different color) in a single plot command without having to manipulate the data. Here's the simplest way I could think of - create a series (column) for each point).

x3 = [x2;x2]

x3 = 2×5

2 3 4 5 6 2 3 4 5 6

y2 = [y;y]

y2 = 2×5

1 2 3 4 5 1 2 3 4 5

figure
plot(x3,y2,'*')

[Guilherme de Melo]

Hi,

Ok, here is example with values in variables. See that I included a high error and it should be plotted in red color from the colorbar:

back_azm_ev2_error=[5, 0.1, 0.5, 0.03, 34];
time_vector_ev2=[600, 605, 610, 615, 620];
baz_lines_ev2=[ 92, 90, 91, 95, 92];

total_elements_baz_lines_ev2=numel(back_azm_ev2_error);
colororder(jet(total_elements_baz_lines_ev2)); %preparing palete color
plot(time_vector_ev2,baz_lines_ev2,'*')
colormap jet
t = colorbar('Ticks',[0 10]);
set(get(t,'ylabel'),'String','Error, degree','fontsize',14);
caxis([0 10]);

[Cris LaPierre]

You create a variable that contains a measure of error, but you don't plot it. How is MATLAB to know that information?

You will need to pass that information into your figure if you want it to be used, and you have to do so in a may that MATLAB interprets as marker color.

I think the simplest way is to use scatter, which allows you to pass in a vector of numbers that is interpretted as colormap indices using the following syntax:

• scatter(x,y,sz,c)

back_azm_ev2_error=[5, 0.1, 0.5, 0.03, 34];
time_vector_ev2=[600, 605, 610, 615, 620];
baz_lines_ev2=[ 92, 90, 91, 95, 92];

total_elements_baz_lines_ev2=numel(back_azm_ev2_error);
colororder(jet(total_elements_baz_lines_ev2)); %preparing palete color
% Note: it is also possible to adjust marker size based on error, but not
% required. Use the default size by passing in an empty array, []
scatter(time_vector_ev2,baz_lines_ev2,[],back_azm_ev2_error,'*')
colormap jet
t = colorbar('Ticks',[0 10]);
set(get(t,'ylabel'),'String','Error, degree','fontsize',14);
caxis([0 10]);

[Guilherme de Melo]

Hi. It worked well. Thank you. I did not know about the scatter. I was thinking that it would be plotted automatically after define the color order, similar to we have done to the lines on map.

However, there is a "problem?" that after use the scatter function, some vertical lines (e.g. 90-95 of y-axis in 605 x-axis ) provide by xlim are not being plotted. They were working normal after "plot" function.

[Cris LaPierre]

If you want to use colororder to set your marker color, see my previous reply. You may also find the description on the Y-input on the plot documenation page helpful.

For your new question, there was no use of xlim in the code you shared. Please provide example code.