%
% Plots the tracks and intensities of events from best track data sets.
% Tracks plotted using Google Earth.
% The latter must be available in matlab binary format. 
% Copyright 2009 WindRiskTech L.L.C.
% Last modified October, 2018
%
%------------------------------------------------------------------
units='knots';      %  Velocities in 'knots' or 'ms' (meters per second)
twidthge=8;         %  Width of track
arrowangle=45;      %  Angle of arrowheads
arrowlength=1;      %  Length of arrowheads
dotsize=1;          %  Size of dots at date positions
dotcolor=[1 0 0];   %  Color of date dots
textoffset=4;       %  Latitude and Longitude offset of text (degrees)
getextsize=1;       %  Scale factor of text size
%------------------------------------------------------------------       
try
    kbestplot=kml('getempfile');
catch err
    errordlg('The free KML toolbox must be installed into MATLAB to use this script')
    return
end      
%
params    %  Load parameters
%
cspline='y';
%
if longmin < 0
   longmin=longmin+360;
end
if longmax < 0
   longmax=longmax+360;
end
%
% Specify storm basin (at, ep, wp, io, sh, or bm)
%
	basin=input('Input basin (al, ep, wp, io, sh, or bm)','s');
%
basin2=basin;
if strcmp(basin,'al')
   basin2='at';
end   
if strcmp(basin,'bm')
    basin2='bom';
end    
file=strcat(trkdir,'/',basin2,'tracks');
atemp=units;
clear x* y* v* u* date offset
units=atemp;
xmin=longmin;
xmax=longmax;
ymin=latmin;
ymax=latmax;
load(file)  
[iq,mq]=size(daym);
year=input('Year of storm');
clear temp
temp=input('Storm number or name','s');
clear ns
ns=str2num(temp); %#ok<ST2NM>
if isempty(ns)
    for i=1:iq
        if strcmpi(temp,strtrim(stname(i,:))) && yeari(i) == year
            ns=nsi(i);
        end
    end
    if isempty(ns)
        disp('Storm not found, quitting')
        return
    end
end      
ystring=num2str(year);
ystringp=num2str(year+1);
ystringm=num2str(year-1);
nstring=num2str(ns);
months={'January','February','March','April','May','June','July','August','September','October','November','December'};
colormap('default')
%colormap(trackcolormap)
colormap=jet(64);
vcolor=colormap;
[qcolor,~]=size(vcolor);
n=0;
m=0;
iswitch=0;
for i=1:iq
   if year == yeari(i) && ns == nsi(i)
      clear stormname
      if min(size(stname)) == 1
          stormname=stname(i);
      else    
        stormname=strtrim(stname(i,:)); 
      end  
      x=zeros(1,mq);
      y=x;
      v=x;
      xd=x;
      yd=x;
      day=x;
      month=cell(mq);
      date=cell(mq);
      for j=1:mq
         if latm(i,j)~=0
            iswitch=1;
            n=n+1;
            x(n)=longm(i,j);
            if strcmp(basin,'at') && x(n) < 100
               x(n)=x(n)+360;
            end   
            y(n)=latm(i,j);
            v(n)=vsm(i,j);
            if hourm(i,j)<2
               m=m+1;
               xd(m)=x(n);
               yd(m)=y(n);
               day(m)=daym(i,j);
               month(m)=cellstr(num2str(monthm(i,j)));
               daystring=num2str(day(m));
               date(m)=cellstr(daystring);
            end      
         end
      end
      x(n+1:end)=[];
      y(n+1:end)=[];
      v(n+1:end)=[];
      xd(m+1:end)=[];
      yd(m+1:end)=[];
      day(m+1:end)=[];
      month(m+1:end)=[];
      date(m+1:end)=[];
   end
end
if iswitch == 0
   message='No data found for this storm, quitting';
   disp(message)
	return   
end
%
if strcmp(cspline,'y')
   b=max(size(x));
   bx=1:b;
   tint=floor(6/trackres);
   tintin=1/tint;
   bxm=1:tintin:b;
   x=spline(bx,x,bxm);
   y=spline(bx,y,bxm);
   v=spline(bx,v,bxm);
end   
a=size(x);
%
x=wrapTo180(x);
xd=wrapTo180(xd);
%
b=a(1,2);
if colorcoded > 0
   for i=2:b
      xt(1)=x(i-1);
      yt(1)=y(i-1);
      xt(2)=x(i);
      yt(2)=y(i);
      vta=v(i-1);
      if colorcoded == 1
         if vta < 33
            vtcolor=vcolor(1,:);
         elseif vta >= 33 && vta < 45
            vtcolor=vcolor(13,:);
         elseif vta >= 45 && vta < 54.4
            vtcolor=vcolor(30,:);
         elseif vta >= 54.4 && vta < 62.4
            vtcolor=vcolor(39,:);
         elseif vta >= 62.4 && vta < 69.5
            vtcolor=vcolor(49,:);
         else
            vtcolor=vcolor(61,:);
         end
      else   
         vtr=round(vta)-15;
         if vtr > 64
            vtr=64;
         end
         if vtr < 1
            vtr=1;
         end   
         vtcolor=vcolor(vtr,:);
      end   
      sar=vtcolor;
      scolor=(strcat('FF',rgbconv(fliplr(vtcolor))));
      kbestplot.plot(xt,yt,'linecolor',scolor,'linewidth',twidthge);
   end   
else
   scolor=(strcat('FF',rgbconv(fliplr(sar))));
   kbestplot.plot(x,y,lc,'linecolor',scolor,'linewidth',twidthge)
   xt(1)=x(end-1);
   xt(2)=x(end);
   yt(1)=y(end-1);
   yt(2)=y(end);
end
%
if strcmp(cplot,'y')
   sd=strcat('FF',rgbconv(fliplr(dotcolor)));
   for jj=1:max(size(xd))
        kbestplot.scatter(xd(jj),yd(jj),'iconcolor',sd,'iconscale',dotsize); 
   end     
end 
%
bd=size(xd);
bd=bd(1,2);
if dplot == 'y'
    if bd > 1   
       ul=zeros(1,bd);
       vl=ul;
       for kl=2:bd-1
          ul(kl)=0.5.*(xd(kl+1)-xd(kl-1));
          vl(kl)=0.5.*(yd(kl+1)-yd(kl-1));
       end
       ul(1)=4.*(x(2)-x(1));
       vl(1)=4.*(y(2)-y(1));
       ul(bd)=4.*(x(b)-x(b-1));
       vl(bd)=4.*(y(b)-y(b-1));
    else
       ul(1)=x(2)-x(1);
       vl(1)=y(2)-y(1);
    end   
    offset=zeros(1,bd)-0.2;
    for kl=1:bd
       if day(kl)>9
          offset(kl)=-0.5;
       end
    end   
    vabs=0.001+sqrt(ul.^2+vl.^2);
    mfactorx=(xmax-xmin)/30;
    mfactory=(ymax-ymin)/30;
    dx=mfactorx*textdfactor*((vl./vabs)+offset);
    dy=mfactory*textdfactor*(-ul./vabs);
    xdb=xd+dx;
    ydb=yd+dy;
    %scolor=(strcat('FF',rgbconv(fliplr(textcolor))));
    for jj=1:max(size(xdb))
        ud=strcat(month(jj),'/',date(jj));
        kbestplot.text(xdb(jj),ydb(jj),0,ud,'labelscale',getextsize);
    end    
end  
%
xquiv=xt(2);
yquiv=yt(2);
uquiv=(xt(2)-xt(1));
vquiv=(yt(2)-yt(1));    
di=1./sqrt(uquiv^2+vquiv^2)+1e-8;
ct=uquiv.*di;
st=vquiv.*di;
angle=pi*(180-arrowangle)/180;
san=sin(angle);
can=cos(angle);
dx1=arrowlength*(ct*can+st*san);
dx2=arrowlength*(ct*can-st*san);
dy1=arrowlength*(st*can-ct*san);
dy2=arrowlength*(st*can+ct*san);
xa=[xt(2) xt(2)+dx1];
ya=[yt(2) yt(2)+dy1];
xb=[xt(2) xt(2)+dx2];
yb=[yt(2) yt(2)+dy2];
kbestplot.plot(xa,ya,'linecolor',scolor,'linewidth',twidthge); 
kbestplot.plot(xb,yb,'linecolor',scolor,'linewidth',twidthge);  
%
if colorcoded > 0 && strcmp(cbar,'y')
   if colorcoded == 1
      ucolor=[vcolor(1,:);vcolor(13,:);vcolor(30,:);vcolor(39,:);vcolor(49,:);vcolor(61,:)];
      aa=0:6;
      bb=1/6;
      tstring='Category';
   else
      ucolor=vcolor;
      if strcmp(units,'ms')
          aa=20:10:80;
          bb=1/6;
          tstring='m/s';
      else
          aa=40:20:160;
          bb=1/6;
          tstring='knots';
      end    
   end   
   kbestplot.colorbar([0,max(aa)],ucolor);
end
%
kbestplot.run
delete('Colorbar.png')