%
%  This script creates maps of the difference between simulated and
%  observed translation velocity vectors. It also calculates and displays
%  the best-track-cell-count-weigthed means of the zonal and meridional
%  translation velocity differences between simulated and observed tracks
%  as well as the weighted root mean square translation velocity
%  difference.The quantity mres specifies the resolution of the maps. This 
%  routine requires installation of the "m_map" routines.
%
%  April, 2017
%  Copyright 2017 WindRiskTech L.L.C. 
%  Last modified June, 2022 to include month filter
%-------------------------------------------------------------------------
months=1:12;
%-------------------------------------------------------------------------
%
params    %  Load parameters
%
% Determine if MATLAB Mapping Toolbox is present
%
v=ver;
maptoolpresent=any(strcmp('Mapping Toolbox', {v.Name}));
clear v
if strcmp(MapOverride,'y')
    maptoolpresent=0;
end   
%
if exist('shearstore','var') == 0 || exist('yearstore','var') == 0
    shape='circ';
    load temp
end
load sorted
clf('reset')
clear ncells zx2 zy2 ncsum
%
kmfac=60*1.852*0.5;
kmfacb=kmfac;
xbound=0;
if strcmp(bas,'AL')
    xbound=260;
end    
%
[n,m]=size(vnet);
%
[~,jmax]=min(vnet,[],2);
jmax=jmax-1;
if exist('yearstore','var') == 0
    yearstore=2000+zeros(1,n);
    freqyear=freq;
end
%
%  Implement velocity, latitude, longitude, and month bounds on calculations
%
latmask=min(ceil(abs(latstore)),1);
vmask=min(ceil(max(vnet-vmaxdisp,0)),1);
monthmask=ismember(monthstore,months);
netmask=latmask.*vmask.*monthmask;
if strcmp(tauto, 'n')
    latmask=latmask.*min(ceil(max(abs(latstore)-tlatmin,0)),1).*min(ceil(max(tlatmax-abs(latstore),0)),1);
    longmask=min(ceil(max(longstore-tlongmin,0)),1).*min(ceil(max(tlongmax-longstore,0)),1);
    netmask=latmask.*longmask.*vmask.*monthmask;
end    
%
if strcmp(bas,'MT') || max(latstore(:,1)) >= 55   %  For Mediterranean basin, or polar lows rectify longitudes
    for i=1:n
        for j=1:m
            if longstore(i,j) >= 200
                longstore(i,j)=longstore(i,j)-360.001; %#ok<SAGROW>
            end
        end    
    end
end      
%
pi=acos(-1);
pifac=pi/180;
fac=2*3600*(1852/3600)^3; % Factor to convert power dissipation rate to m^3/s^2
fac=fac/mres^2;
%
if strcmp(bas,'GB')
   projection=gproject;
end   
%
if strcmp(mapmode,'auto')
    xmin=min(nonzeros(longstore));
    xmax=max(nonzeros(longstore));
    ymin=min(nonzeros(latstore));
    ymax=max(nonzeros(latstore));
    xmin=mres*floor(xmin/mres);
    xmax=mres*ceil(xmax/mres);
    ymin=mres*floor(ymin/mres);
    ymax=mres*ceil(ymax/mres);
else
    xmin=longmin;
    xmax=longmax;
    ymin=latmin;
    ymax=latmax;
end 
if strcmp(bas,'GB')
    xmin=0;
    xmax=360;
    ymin=-87;
    ymax=87;
end 
xmin2=xmin-mres;
xmax2=xmax+mres;
ymin2=ymin-mres;
ymax2=ymax+mres;
%
y=ymin2:mres:ymax2;
x=xmin2:mres:xmax2; 
%
nx=max(size(x));
ny=max(size(y));
%
yearmin=min(yearstore);
yearmax=max(yearstore);
qy=yearmax-yearmin+1;
year=yearmin:yearmax;
longstore=longstore+rand*1e-6; % Added September 2020 to avoid actual zeros in longitude
%
dum=nonzeros(longstore.*netmask);
longr=mres*floor(dum/mres);
longr=longr+mod(x(1),mres);
clear dum
%
dum=nonzeros(netmask.*latstore);
latr=mres*floor(dum/mres);
latr=latr+mod(y(1),mres);
%
delx=zeros(n,m);
dely=delx;
%
% Correct for tracks crossing the Greenwich Meridian
%
for i=1:n
    u=longstore(i,2:jmax(i))-longstore(i,1:jmax(i)-1);
    sdelx=sign(u).*sign(180-abs(u));
    ua=min(abs(u),abs(abs(u)-360));
    delx(i,2:jmax(i))=cos(pifac.*latstore(i,2:jmax(i))).*ua.*sdelx;
    dely(i,2:jmax(i))=latstore(i,2:jmax(i))-latstore(i,1:jmax(i)-1);
end
clear u ua sdelx
%
delxr=nonzeros(netmask.*(delx+0.001*pi))-0.001*pi;
delyr=nonzeros(netmask.*(dely+0.001*pi))-0.001*pi;
dum=repmat(max(yearstore',1),[1,m]); 
yr=nonzeros(double(dum).*netmask);
if max(size(freqyear)) > 1
    atemp=ismember((min(yr):max(yr)),yr);
    atemp=cast(atemp,'like',yr);
    yearset=nonzeros(atemp.*(min(yr):max(yr)));
    [~,nyy]=ismember(yr,yearset);
    freqmask=freqyear(nyy);
else
    freqmask=freq;
end    
clear dum atemp
%
[latl,ay]=ismember(latr,y);
[longl,ax]=ismember(longr,x);
[yearl,an]=ismember(yr,year);
ax=max(ax,1);
ay=max(ay,1);
an=max(an,1);
%
zx=zeros(nx,ny,qy);
zy=zeros(nx,ny,qy);
ncells=zeros(nx,ny,qy);
N=max(size(latl));
delxnet=kmfac.*delxr.*freqmask'.*latl.*longl.*yearl./freq;
delynet=kmfac.*delyr.*freqmask'.*latl.*longl.*yearl./freq;
for i=1:N
    zx(ax(i),ay(i),an(i))=zx(ax(i),ay(i),an(i))+delxnet(i);
    zy(ax(i),ay(i),an(i))=zy(ax(i),ay(i),an(i))+delynet(i);
    ncells(ax(i),ay(i),an(i))=ncells(ax(i),ay(i),an(i))+1;
end
% 
zx2(:,:)=mean(zx./(ncells+1e-6),3);
zy2(:,:)=mean(zy./(ncells+1e-6),3);
ncsum(:,:)=sum(ncells,3);
%
%  Now process best tracks
%
bestproc   
%
yearbeg=min(bestyears);
yearend=max(bestyears);
yearsize=max(size(bestyears));
%
[nb,mb]=size(vbest);
%
%  Implement velocity, latitude, and longitude bounds on calculations
%
latmaskb=min(ceil(abs(latbest)),1);
vmaskb=min(ceil(max(vbest-vmaxdisp,0)),1);
dum=datevec(datebest);
monthbest=dum(:,2);
monthbest=reshape(monthbest,[nb,mb]);
monthmaskb=ismember(monthbest,months);
netmaskb=latmaskb.*vmaskb.*monthmaskb;
if strcmp(tauto, 'n')
    latmaskb=latmaskb.*min(ceil(max(abs(latbest)-tlatmin,0)),1).*min(ceil(max(tlatmax-abs(latbest),0)),1);
    longmaskb=min(ceil(max(longbest-tlongmin,0)),1).*min(ceil(max(tlongmax-longbest,0)),1);
    netmaskb=latmaskb.*longmaskb.*vmaskb.*monthmaskb;
end    
%
[~,jmaxb]=min(vbest,[],2);
jmaxb=jmaxb-1;
%
dum=nonzeros(longbest.*netmaskb);
longrb=mres*floor(dum/mres);
longrb=longrb+mod(x(1),mres);
clear dum
%
dum=nonzeros(netmaskb.*latbest);
latrb=mres*floor(dum/mres);
latrb=latrb+mod(y(1),mres);
%
[latlb,ayb]=ismember(latrb,y);
[longlb,axb]=ismember(longrb,x);
axb=max(axb,1);
ayb=max(ayb,1);
%
delxb=zeros(nb,mb);
delyb=delxb;
%
% Correct for tracks crossing the Greenwich Meridian
%
for i=1:nb
    u=longbest(i,2:jmaxb(i))-longbest(i,1:jmaxb(i)-1);
    sdelx=sign(u).*sign(180-abs(u));
    ua=min(abs(u),abs(abs(u)-360));
    delxb(i,2:jmaxb(i))=cos(pifac.*latbest(i,2:jmaxb(i))).*ua.*sdelx;
    delyb(i,2:jmaxb(i))=latbest(i,2:jmaxb(i))-latbest(i,1:jmaxb(i)-1);
end
clear u ua sdelx
%for i=1:nb,
%    delxb(i,2:jmaxb(i))=cos(pifac.*latbest(i,2:jmaxb(i))).*(longbest(i,2:jmaxb(i))-longbest(i,1:jmaxb(i)-1));
%    delyb(i,2:jmaxb(i))=latbest(i,2:jmaxb(i))-latbest(i,1:jmaxb(i)-1);
%end
delxrb=nonzeros(netmaskb.*(delxb+0.01*pi))-0.01*pi;
delyrb=nonzeros(netmaskb.*(delyb+0.01*pi))-0.01*pi;
%
zxb=zeros(nx,ny);
zyb=zeros(nx,ny);
ncb=zeros(nx,ny);
Nb=max(size(latlb));
delxnetb=kmfacb.*delxrb.*latlb.*longlb;
delynetb=kmfacb.*delyrb.*latlb.*longlb;
for i=1:Nb
    zxb(axb(i),ayb(i))=zxb(axb(i),ayb(i))+delxnetb(i);
    zyb(axb(i),ayb(i))=zyb(axb(i),ayb(i))+delynetb(i);
    ncb(axb(i),ayb(i))=ncb(axb(i),ayb(i))+1;
end
% 
zzxb=(zxb./(ncb+1e-6));
zzyb=(zyb./(ncb+1e-6));
%
% Calculate new map limits
% 
imin=1;  
imax=ny;
kmin=1;
kmax=nx;
if strcmp(mapmode,'auto')
    xmin1=xmax;
    xmax1=xmin;
    ymin1=ymax;
    ymax1=ymin;
    kmin=nx;
    kmax=1;
    imin=ny;
    imax=1;
    zcrit=0.1*max(abs(zx2));
    for j=2:ny-1
        for k=2:nx-1 
            if abs(zx2(k,j)) > zcrit & x(k) > xbound %#ok<AND2>
                xmin1=min(xmin1,x(k));
                xmax1=max(xmax1,x(k));
                ymin1=min(ymin1,y(j));
                ymax1=max(ymax1,y(j));
                kmin=min(kmin,k);
                kmax=max(kmax,k);
                imin=min(imin,j);
                imax=max(imax,j);
            end    
        end
    end    
    xmin=xmin1;
    ymin=ymin1;
    xmax=xmax1;
    ymax=ymax1;
end
%
% Calculate biases and RMS difference between best track and synthetic data
%
usyn=zx2(kmin:kmax,imin:imax);
vsyn=zy2(kmin:kmax,imin:imax);
ubest=zzxb(kmin:kmax,imin:imax);
vbest=zzyb(kmin:kmax,imin:imax);
udif=usyn-ubest;
vdif=vsyn-vbest;
ncmap=ncb(kmin:kmax,imin:imax);
ncmap(ncmap<mintracks)=0; % Assign zero weight where best tracks are deficient
kdif=kmax-kmin+1;
idif=imax-imin+1;
zxdif=reshape(udif,[1,kdif*idif]);
zydif=reshape(vdif,[1,kdif*idif]);
ncweight=reshape(ncmap,[1,kdif*idif]);
rmsdif=sqrt(sum(ncweight.*(zxdif.^2+zydif.^2))/sum(ncweight));  % RMS weighted by best track counts
absdif=sum(ncweight.*sqrt(zxdif.^2+zydif.^2))/sum(ncweight);    % Mean absolute value weighted by best track counts
xbias=sum(ncweight.*zxdif)/sum(ncweight);  % Zonal bias weighted by best track counts
ybias=sum(ncweight.*zydif)/sum(ncweight);  % Meridional bias weighted by best track counts
xbias=0.01*round(100*xbias);  % Round to three significant figures
ybias=0.01*round(100*ybias);
rmsdif=0.01*round(100*rmsdif);
%
label1=['X-bias = ',num2str(xbias),' km/hr'];
label2=['Y-bias = ',num2str(ybias),' km/hr'];
label3=['RMS difference = ',num2str(rmsdif),' km/hr'];
mtit='Difference Stats';
fm=msgbox({label1;label2;label3},mtit);
set(fm,'resize','on')
%
%  Calculate and display difference vectors
%
clf('reset')
%
% Set up map
%
%  Condition translation vectors and maximum magnitudes
%
udif(ncmap<mintracks)=NaN;  % Do not display differences where there are insufficient historical tracks
vdif(ncmap<mintracks)=NaN;
usyn(ncsum(kmin:kmax,imin:imax)<mintracks)=NaN;
vsyn(ncsum(kmin:kmax,imin:imax)<mintracks)=NaN;
ubest(ncmap<mintracks)=NaN;
vbest(ncmap<mintracks)=NaN;
%
difmax=round(nanmax(nanmax(sqrt(udif.^2+vdif.^2))));
synmax=round(nanmax(nanmax(sqrt(usyn.^2+vsyn.^2))));
bestmax=round(nanmax(nanmax(sqrt(ubest.^2+vbest.^2))));
%
[x2,y2]=meshgrid(x(kmin:kmax),y(imin:imax));
%
if strcmp(bas,'GB')
    projection=gproject;
end  
%
if maptoolpresent
    %gh=figure('Visible','off');
    %
    axesm('MapProjection',projection,'MapLatLimit',[ymin ymax],'MapLonLimit',[xmin xmax], ...
         'frame','on','ffacecolor',oceancolor,'fontname',axisfont,'fontsize',axisfontsize, ...
         'fontweight',axisfontweight,'labelrotation','on')
    delx=floor((xmax-xmin)/5);
    dely=floor((ymax-ymin)/5);
    gridm('mlinelocation',delx,'plinelocation',dely,'mlabellocation',delx,'plabellocation',dely,'mlabelparallel','south')
    if strcmp(gridline,'none')
        gridm('off')
    end
    plabel('fontweight','bold'); mlabel('fontweight','bold');
    axis off
    %R = georasterref('RasterSize', size(z(imin:imax,kmin:kmax)), ...
    %   'Latlim', [ymin ymax], 'Lonlim', [xmin xmax]);
    R = georasterref('RasterSize', [imax-imin+1, kmax-kmin+1],'Latlim', [ymin ymax], 'Lonlim', [xmin xmax]);
    %
    geoback(xmin, xmax, ymin, ymax)
    axis off
    mlabel('south')
    %
    % Map title
    %
    if isempty(storm)
        datea(1)=upper(datea(1));
        tit1=datea;
    else
        tit1=[storm(1:2) ' ' storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
        ' ' datea(9:10) ' GMT'];
    end 
    %
    tightmap
    im=2;
    while im > 1
        im=menu('Choose Plot','Stop','Synthetic','Historical','Both','Difference');
        delete(fm)
        if im > 1
            try
                delete(hb)
            catch    
            end
            try
                delete(h2b)
            catch    
            end    
            try 
                delete(an)
            catch
            end            
        end
        if im == 1
            break
        elseif im == 2
            hb=quiverm(y2,x2,vsyn',usyn',quiv_color,2); 
            tita='Mean Synthetic Track Vectors';
            tittot={tita;tit1};
            title(tittot,'fontsize',12,'fontweight','bold','interpreter','none')
            %
            % Display magnitudes in text box
            %
            dim = [0.64 0 0.3 0.1]; % Box position
            str1=['Maximum magnitude = ',num2str(synmax),' km/hr'];
            an=annotation('textbox',dim,'String',str1,'FitBoxToText','on'); 
        elseif im == 3
            h2b=quiverm(y2,x2,vbest',ubest','b',2);
            tita='Mean Historical Track Vectors';
            tittot={tita;tit1};
            title(tittot,'fontsize',12,'fontweight','bold','interpreter','none')
            %
            % Display magnitudes in text box
            %
            dim = [0.64 0 0.3 0.1]; % Box position
            str1=['Maximum magnitude = ',num2str(bestmax),' km/hr'];
            an=annotation('textbox',dim,'String',str1,'FitBoxToText','on');             
        elseif im == 4
            hb=quiverm(y2,x2,vsyn',usyn',quiv_color,2);
            h2b=quiverm(y2,x2,vbest',ubest','b',2);
            tita='Mean Synthetic and Historical Track Vectors';
            tittot={tita;tit1};
            title(tittot,'fontsize',12,'fontweight','bold','interpreter','none')            
        elseif im == 5   
            hb=quiverm(y2,x2,vdif',udif',quiv_color,2);
            tita='Mean Track Difference Vectors';
            tittot={tita;tit1};
            title(tittot,'fontsize',12,'fontweight','bold','interpreter','none')    
            %
            % Display magnitudes in text box
            %
            dim = [0.6 0.02 0.3 0.1]; % Box position
            str1=['Maximum magnitude = ',num2str(difmax),' km/hr'];
            str2=['RMS difference = ',num2str(0.1*round(10*rmsdif)),' km/hr'];
            str={str1;str2};
            an=annotation('textbox',dim,'String',str,'FitBoxToText','on');             
        end  
        tightmap
        %
        % Add button to access documentation of script
        %
        if strcmp(Docdisp,'y')
            try
                delete(hButton)
            catch
            end    
            tempname=mfilename;
            mname = strcat('graphics_docs/html/',tempname,'_doc.html');
            comstr=['open ',mname,';delete(hButton)'];
            labelStr = '<html><center><a href="">Doc';
            hButton = uicontrol('string',labelStr,'pos',[10,10,50,20],'callback',comstr); 
        end                          
    end    
    hold off
else    
    m_proj(projection,'long',[xmin xmax],'lat',[ymin ymax]);
    m_gridm('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight, ...
        'linestyle',gridline,'backgroundcolor',oceancolor)
    if strcmp(bas,'MT')
         m_gshhs_l('patch',landcolor,'edgecolor','none');
    else     
         m_coast('patch',landcolor,'edgecolor','none');
    end   
    if strcmp(pstates,'y')
        m_states(xmin, xmax, ymin, ymax,'color',stcolor)
    end 
     %
    % Map title
    %
    if isempty(storm)
        datea(1)=upper(datea(1));
        tit1=datea;
    else
        tit1=[storm(1:2) ' ' storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
        ' ' datea(9:10) ' GMT'];
    end  
    %
    hold on
    im=2;
    while im > 1
        im=menu('Choose Plot','Stop','Synthetic','Historical','Both','Difference');
        delete(fm)
        if im > 1
            try
                delete(hb)
            catch    
            end
            try
                delete(h2b)
            catch    
            end    
            try 
                delete(an)
            catch
            end
        end
        if im == 1
            break
        elseif im == 2
            hb=m_quiver(x2,y2,usyn',vsyn',2,'color',quiv_color);
            tita='Mean Synthetic Track Vectors';
            tittot={tita;tit1};
            title(tittot,'fontsize',12,'fontweight','bold','interpreter','none')
            %
            % Display magnitudes in text box
            %
            dim = [0.64 0 0.3 0.1]; % Box position
            str1=['Maximum magnitude = ',num2str(synmax),' km/hr'];
            an=annotation('textbox',dim,'String',str1,'FitBoxToText','on'); 
        elseif im == 3
            h2b=m_quiver(x2,y2,ubest',vbest',2,'color','b');
            tita='Mean Historical Track Vectors';
            tittot={tita;tit1};
            title(tittot,'fontsize',12,'fontweight','bold','interpreter','none')
            %
            % Display magnitudes in text box
            %
            dim = [0.64 0 0.3 0.1]; % Box position
            str1=['Maximum magnitude = ',num2str(bestmax),' km/hr'];
            an=annotation('textbox',dim,'String',str1,'FitBoxToText','on'); 
        elseif im == 4
            hb=m_quiver(x2,y2,usyn'./10,vsyn'./10,0,'color',quiv_color);
            h2b=m_quiver(x2,y2,ubest'./10,vbest'./10,0,'color','b');
            tita='Mean Synthetic and Historical Track Vectors';
            tittot={tita;tit1};
            title(tittot,'fontsize',12,'fontweight','bold','interpreter','none')
        elseif im == 5    
            hb=m_quiver(x2,y2,udif',vdif',2,'color',quiv_color);
            tita='Mean Track Difference Vectors';
            tittot={tita;tit1};
            title(tittot,'fontsize',12,'fontweight','bold','interpreter','none')   
            %
            % Display magnitudes in text box
            %
            dim = [0.6 0.02 0.3 0.1]; % Box position
            str1=['Maximum magnitude = ',num2str(difmax),' km/hr'];
            str2=['RMS difference = ',num2str(0.1*round(10*rmsdif)),' km/hr'];
            str={str1;str2};
            an=annotation('textbox',dim,'String',str,'FitBoxToText','on');             
        end  
        %
        % Add button to access documentation of script
        %
        if strcmp(Docdisp,'y')
            try
                delete(hButton)
            catch
            end    
            tempname=mfilename;
            mname = strcat('graphics_docs/html/',tempname,'_doc.html');
            comstr=['open ',mname,';delete(hButton)'];
            labelStr = '<html><center><a href="">Doc';
            hButton = uicontrol('string',labelStr,'pos',[10,10,50,20],'callback',comstr); 
        end                          
    end    
    hold off
end
%
if strcmp(bas,'MT') %  For Mediterranean basin, reset longitudes
    load temp longstore
end  
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 