clc
clear all
close all

% GPS positions (latitude,longitude) of some markers
data = [
    -22.976730, - 43.195080 ;
     55.756950,   37.614975 ;
     33.605381, -  7.631940 ;
     35.670479,  139.740921 ;
     51.506325, -  0.127144 ;
     40.714550, - 74.007124 ;
    -33.869629,  151.206955 ;
    -26.204944,   28.040035 ;
     37.777125, -122.419644 ;
     30.083740,   31.255360 ;
      6.439180,    3.423480
];
labels = {
    'Rio de Janeiro'
    'Moscow'
    'Casablanca'
    'Tokyo'
    'London'
    'New York'
    'Sydney'
    'Johannesburg'
    'San Francisco'
    'Cairo'
    'Lagos'
};

storm_name = 'Idai';
start_year = 2019;

track_file = ['/discover/nobackup/jkolassa/SMAP_TC_project/track_files/' storm_name '_track.ascii'];
track_data = dlmread(track_file);
track_datetime = datetime(start_year, track_data(:,1), track_data(:,2), track_data(:,3), 0, 0);
track_datenum = datenum(start_year, track_data(:,1), track_data(:,2), track_data(:,3), 0, 0);
track_lat = track_data(:,5);
track_lon = track_data(:,4);

data = [track_lat track_lon];

track_datenum = round(track_datenum);
unique_datenum = unique(track_datenum);
color_vec = NaN*ones(length(track_datenum),3);

cmap = colormap(jet(length(unique_datenum)));

for i = 1:length(unique_datenum)

    date_index = find(track_datenum==unique_datenum(i));
    for j = 1:length(date_index)
        color_vec(date_index(j),:) = cmap(i,:);
    end
end

% world map in Mercator projection
fname = 'snapshot-2022-11-29T00_00_00Z_small.jpeg';
img = imread(fname);
[imgH,imgW,~] = size(img);

% Mercator projection
[x,y] = mercatorProjection(data(:,2), data(:,1), imgW, imgH);

% plot markers on map
imshow(img, 'InitialMag',100, 'Border','tight'), hold on
for i = 1:length(unique_datenum)
date_index = find(track_datenum==unique_datenum(i));
plot(x(date_index),y(date_index), 'o', 'MarkerSize',5, 'MarkerFaceColor', cmap(i,:), 'MarkerEdgeColor', cmap(i,:), 'LineWidth',0.5)
end
%text(x, y, track_datetime, 'Color','w', 'VerticalAlign','bottom', 'HorizontalAlign','right')
hold off
saveas(gcf,'test_map_2.png')

function [x,y] = mercatorProjection(lon, lat, width, height)
    x = mod((lon-24)*width/28, width) ;
%    y = height/2 - log(tan((lat+90)*pi/360))*width/(2*pi);
   y = mod((90-100-lat)*height/14, height) ;
end
