clear
prtflag   = 0; % flag
figdir    = '/Users/hildeoliver/Documents/rgb19/figs';
matdir    = '/Users/hildeoliver/Documents/rgb19/scripts/sst/mat/';
%files     = dir(matdir);
%files.name = ['udelsst_mab_20190515_2135.mat'];
%files.name = ['udelsst_mab_20190515_2348.mat'];
%files.name = ['udelsst_mab_20190516_0128.mat'];
%files.name = ['udelsst_mab_20190516_2124.mat'];
%files.name = ['udelsst_mab_20190516_2337.mat'];
%files.name = ['udelsst_mab_20190517_0947.mat'];
%files.name = ['udelsst_mab_20190517_1338.mat'];
%files.name = ['udelsst_mab_20190518_1327.mat'];
%files.name = ['udelsst_mab_20190519_0052.mat'];
files.name = ['udelsst_mab_20190520_0041.mat'];

% bird data
birdfile = '/Users/hildeoliver/Documents/rgb19/RB1904_birds.csv';
birdtable = readtable(birdfile);
birdid = repmat('o',1, size(birdtable,1));
% whales and dolphins
wod = strcmp(birdtable.SpeciesCode, 'CODO') | strcmp(birdtable.SpeciesCode, 'UNWH') ...
     | strcmp(birdtable.SpeciesCode, 'SPWH') | strcmp(birdtable.SpeciesCode, 'UNGD') ...
     | strcmp(birdtable.SpeciesCode, 'UNDO') | strcmp(birdtable.SpeciesCode, 'SEWH');
% gulls
gull = strcmp(birdtable.SpeciesCode, 'HERG') | strcmp(birdtable.SpeciesCode, 'GBBG') ...
     | strcmp(birdtable.SpeciesCode, 'UNGL') | strcmp(birdtable.SpeciesCode, 'LGGB') ...
     | strcmp(birdtable.SpeciesCode, 'LAGU');     
% shearwaters and petrels
sop = strcmp(birdtable.SpeciesCode, 'WISP') | strcmp(birdtable.SpeciesCode, 'MASH') ...
     | strcmp(birdtable.SpeciesCode, 'SOSH') | strcmp(birdtable.SpeciesCode, 'AUSH') ...
     | strcmp(birdtable.SpeciesCode, 'COSH');     
% gannets
gan = strcmp(birdtable.SpeciesCode, 'NOGA');    
% terns
tern = strcmp(birdtable.SpeciesCode, 'COTE') | strcmp(birdtable.SpeciesCode, 'UNST');     
% fulmars
ful = strcmp(birdtable.SpeciesCode, 'NOFU');    
% other
other = ~wod & ~gull & ~sop & ~gan & ~tern & ~ful;

% area for plotting
plonlim = [-71.5 -70.5];
platlim = [39.5 40.5];    % northeast

clrlim = [4 19];

bdata   = load('/Users/hildeoliver/Documents/rgb19/scripts/ngdc2.mat');
blonind = find(bdata.lon>plonlim(1) & bdata.lon < plonlim(2));
blatind = find(bdata.lat>platlim(1) & bdata.lat < platlim(2));
blon    = bdata.lon(blonind);
blat    = bdata.lat(blatind);
bdpth   = bdata.h(blatind(1):blatind(end),blonind(1):blonind(end));

% mooring locations
mrlon  = [-70.88333 -70.88333 -70.775   -70.8833 -70.88333 -70.775 -70.775];
mrlat  = [40.363333 40.226667 40.136667 40.09833 39.94     40.365  39.94];
mrname = ['IS','CI','CN','CO','OS','UI','UO'];
sname  = {'A1','A2','A3','A4','A5','A6','A7','A8','A9','A10','A11','A12','A13','A14','A15','A16','A17','A18'};
zdlon  = [-70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830 -70.830];
zdlat  = [ 40.725  40.660  40.595  40.530  40.465  40.400  40.335  40.270  40.205  40.140  40.075  40.010  39.945  39.880  39.815  39.750  39.685  39.620];
% to read eddy chasing stations
eddypwf = '/Users/hildeoliver/Documents/rgb19/scripts/eddy_waypoints_3.txt';
fid     = fopen(eddypwf);
wpind   = 0;
while ~feof(fid) %Check every line in the file
    wpind = wpind + 1;
    [llstr(wpind,:) count] = fscanf(fid,'%c%c%d    %d %f %c    %d %f %c\n',[1,9]);
end
fclose(fid);
for iw = 1:wpind
    ecsname{iw} = [char(llstr(iw,1)) char(llstr(iw,2)) num2str(llstr(iw,3))];
    latdgr(iw)  = llstr(iw,4);
    latmnt(iw)  = llstr(iw,5);
    latnsl(iw)  = llstr(iw,6);
    londgr(iw)  = llstr(iw,7);
    lonmnt(iw)  = llstr(iw,8);
    lonewl(iw)  = llstr(iw,9);
    ecslat(iw)  = latdgr(iw) + latmnt(iw)./60;
    ecslon(iw)  = londgr(iw) + lonmnt(iw)./60;
    if latnsl(iw) == 83
        ecslat(iw) = -ecslat(iw);
    elseif latnsl(iw) ~= 78
        errmes = ['Latitude letter error in file: ' eddypwf];
        error(errmes);
    end
    if lonewl(iw) == 87
        ecslon(iw) = -ecslon(iw);
    elseif lonewl(iw) ~= 69
        errmes = ['Latitude letter error in file: ' eddypwf];
        error(errmes);
    end
end

for ifile = 1:length(files)
    
    % find the time index to read (note UDel time is generally several
    % minutes earlier than Rutgers
    fname = files(ifile).name;
    if length(fname) > 10 && strcmp(fname(end-2:end),'mat')
        mcsst  = load([matdir fname]);
        lon    = mcsst.lon;
        lat    = mcsst.lat;
        nmcsst = mcsst.mcsst;
        time   = mcsst.time;
                    
        figure(1);clf;
        orient landscape
        colormap('jet');
        
       % axes('pos',[0.51 0.2 0.4 0.75]);
        pcolorjw(lon,lat,nmcsst);
        caxis(clrlim);
        
        hold on;
        [cc,hh] = contour(blon,blat,bdpth,[40:10:100 200 500 1000 2000 3000]);
        set(hh,'color',[.7 .7 .7],'linewidth',.5);
        [cc,hh] = contour(blon,blat,bdpth,[100 200]);
        set(hh,'color',[.7 .7 .7],'linewidth',.85);
        clabel(cc,hh,'color',[0.5 0.5 0.5],'fontsize',10,'labelspacing',350);
        pnc;
        pmc('fill',[0.5 0.5 0.5],[0.5 0.5 0.5]);
        hclb = colorbar;
        set(hclb,'location','southoutside');
        %set(hclb,'pos',[0.52 0.13 0.38 0.02],'ytick',[clrlim(1):2:clrlim(2)],'tickdir','out','fontsize',15);
        xlabel(hclb,'Temperature','fontsize',15);
        %text(plonlim(2)+0.02,platlim(1)-0.2,'(^oC)','fontsize',15);
        xlim(plonlim);ylim(platlim);
        set(gca,'xtick',[plonlim(1):0.25:plonlim(2)],'ytick',[platlim(1):0.25:platlim(2)]);
        set(gca,'layer','top','tickdir','out','fontsize',15,'yticklabel',[]);
        %plot(mrlon,mrlat,'pk','markersize',6,'linewidth',2);
        plot(zdlon,zdlat,'dk','markersize',6,'linewidth',2, 'color',[0.8 0.8 0.8]);
        plot(ecslon,ecslat,'ok','markersize',6,'linewidth',2, 'color',[0.8 0.8 0.8]);
        tttxt   = ['SST on ' datestr(time)];
            hc1 = scatter(birdtable.Long(wod),birdtable.Lat(wod), 150, 'blue','filled','p');            
            hc2 = scatter(birdtable.Long(gull),birdtable.Lat(gull), 75, 'black','o','LineWidth',.9);            
            hc3 = scatter(birdtable.Long(sop),birdtable.Lat(sop), 75, 'black','*','LineWidth',.9);            
            hc4 = scatter(birdtable.Long(gan),birdtable.Lat(gan), 75, 'black','+','LineWidth',.9);            
            hc5 = scatter(birdtable.Long(tern),birdtable.Lat(tern), 75, 'black','o','LineWidth',.9);            
            hc6 = scatter(birdtable.Long(ful),birdtable.Lat(ful), 75, 'black','x','LineWidth',.9);            
            hc7 = scatter(birdtable.Long(other),birdtable.Lat(other), 75, 'black','^','LineWidth',.9);
            legend([hc1 hc2 hc3 hc4 hc5 hc6 hc7], {'dolphins and whales','gulls','shearwaters + petrels',...
                'gannets', 'terns','fulmars','other'},'Location','southwest')

            tttxt   = [tttxt ' + birds'];
            caxis(clrlim);
        htt = title(tttxt);
        set(htt,'fontsize',14,'fontname','Times');

    end
    
end
