clear

% sstime(1) = datenum([2019 04 03 21 28 00]);
% sstime(2) = datenum([2019 04 04 21 16 00]);
% sstime(3) = datenum([2019 04 07 09 16 00]);
% sstime(4) = datenum([2019 04 07 12 27 00]);
% sstime(5) = datenum([2019 04 11 00 01 00]);
% sstime(6) = datenum([2019 04 16 12 01 00]);
% sstime(7) = datenum([2019 04 16 13 05 00]);
% sstime(8) = datenum([2019 04 18 00 20 00]);
% sstime(9) = datenum([2019 04 24 20 45 00]);
% sstime(10) = datenum([2019 04 25 00 38 00]);
% sstime(11) = datenum([2019 04 29 10 02 00]);
% sstime(12) = datenum([2019 04 29 13 55 00]);
% sstime(13) = datenum([2019 05 01 09 38 00]);
% sstime(14) = datenum([2019 05 06 21 46 00]);
% sstime(15) = datenum([2019 05 06 23 18 00]);
% sstime(16) = datenum([2019 05 06 23 58 00]);
% sstime(17) = datenum([2019 05 07 10 10 00]);
% sstime(18) = datenum([2019 05 07 11 34 00]);
% sstime(1) = datenum([2019 05 07 14 02 00]);
% sstime(1) = datenum([2019 05 07 21 35 00]);
% sstime(2) = datenum([2019 05 07 22 53 00]);
% sstime(3) = datenum([2019 05 07 23 47 00]);
% sstime(4) = datenum([2019 05 08 01 27 00]);
% sstime(5) = datenum([2019 05 08 21 23 00]);
% sstime(6) = datenum([2019 05 08 23 35 00]);
% sstime(7) = datenum([2019 05 09 00 08 00]);
% sstime(8) = datenum([2019 05 09 01 15 00]);
% sstime(9) = datenum([2019 05 09 09 46 00]);
% sstime(1) = datenum([2019 05 09 12 24 00]);
% sstime(2) = datenum([2019 05 09 13 39 00]);
% sstime(1) = datenum([2019 05 09 23 42 00]);
% sstime(1) = datenum([2019 05 11 20 48 00]);
% sstime(1) = datenum([2019 05 15 21 42 00]);
% sstime(1) = datenum([2019 05 15 22 52 00]);
% sstime(2) = datenum([2019 05 15 23 54 00]);
% sstime(1) = datenum([2019 05 16 01 35 00]);
% sstime(1) = datenum([2019 05 16 21 31 00]);
% sstime(1) = datenum([2019 05 16 23 42 00]);
% sstime(1) = datenum([2019 05 17 00 07 00]);
% sstime(1) = datenum([2019 05 17 01 23 00]);
% sstime(1) = datenum([2019 05 17 09 47 00]);
%sstime(1) = datenum([2019 05 17 12 23 00]);
%sstime(1) = datenum([2019 05 17 13 45 00]);
%sstime(1) = datenum([2019 05 18 11 57 00]);
%sstime(1) = datenum([2019 05 18 13 34 00]);
%sstime(1) = datenum([2019 05 18 21 07 00]);
%sstime(2) = datenum([2019 05 18 23 17 00]);
sstime(1) = datenum([2019 05 19 00 59 00]);

addadcp   = 0;
readurl   = 1;  % switch to extract data from the UDel server; otherwise read local data

%url       = 'http://tds.maracoos.org/thredds/dodsC/AVHRR.nc';
url       = 'http://basin.ceoe.udel.edu/thredds/dodsC/avhrr_unfiltered_sst.nc';
figdir    = '~/research/cruise/2019_Brown_0512-0525/data/figs/sst/';
savemat   = 1;  % switch to save mat file
matdir    = '~/research/cruise/2019_Brown_0512-0525/data/figs/sst/mat/';

% ADCP data
tranum   = 14;   % the transect number for adding adcp surface current
quvscale = 1;   % quiver vector scale
uvscale  = 0.2; %
adcpfile = '~/research/cruise/2019_Brown_0512-0525/data/adcp/os75bb_uv_blank8mbin4m.mat';
% read transect start time
addpath ~/research/cruise/2019_Brown_0512-0525/data/scripts/
rb1904_transect_startime

lfile   = '~/research/roms/sbfront/data/satellite/lonlat.mat';

% area for extracting data
lonlim  = [-76 -66];  % northeast
latlim  = [36 42];    % northeast
% area for plotting
plonlim = [-71.75 -69.75];
platlim = [39.25 41.1];    % northeast

clrlim = [4 19];

ldata   = load(lfile);
bdata   = load('~/research/roms/bathymetry/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];


if readurl  % read data from UDel server
    udtime = double(ncread(url,'time'))./86400 + datenum([1970 1 1]);
    %udlon  = ncread(url,'lon');
    %udlat  = ncread(url,'lat');
    load('udel_sst_lonlat.mat');
    lonind = find(udlon>lonlim(1) & udlon < lonlim(2));
    latind = find(udlat>latlim(1) & udlat < latlim(2));
    lonlen = length(lonind);
    latlen = length(latind);
    lon    = udlon(lonind);
    lat    = udlat(latind);
else   % read local data
    matfile = dir(matdir);
    for iff = 3:length(matfile)
        matname = [matfile(iff).folder '/' matfile(iff).name];
        tmptime = load(matname,'time');
        udtime(iff-2) = tmptime.time;
    end
    tmplon  = load(matname,'lon');
    tmplat  = load(matname,'lat');
    lon     = tmplon.lon;
    lat     = tmplat.lat;
end

for ifile = 1:length(sstime)
    
    % find the time index to read (note UDel time is generally several
    % minutes earlier than Rutgers
    ttind = find(udtime <= sstime(ifile));
    tind  = ttind(end);
    time  = udtime(tind);
    if abs(time - sstime(ifile)) < 30/(24*60)
        % only save the data when the UDel time is within half a hour
        % window around the Rutgers time
        if readurl
            mcsst = double(ncread(url,'mcsst',[lonind(1) latind(1) tind],[lonlen latlen 1]))';
            if savemat
                matfile = [matdir 'udelsst_mab_' datestr(udtime(tind),'yyyymmdd_HHMM') '.mat'];
                save(matfile,'mcsst','lon','lat','time');
                disp(['Saved file: ' matfile]);
            end
        else
            matname = [matfile(tind+2).folder '/' matfile(tind+2).name];
            tmcsst  = load(matname,'mcsst');
            mcsst   = tmcsst.mcsst;
        end
        plonind = find(lon > plonlim(1) & lon < plonlim(2));
        platind = find(lat > platlim(1) & lat < platlim(2));
        nlon   = lon(plonind);
        nlat   = lat(platind);
        nmcsst = mcsst(platind(1):platind(end),plonind(1):plonind(end));
        
        figure(1);clf;
        orient PORTRAIT
        colormap('jet');
        
        axes('pos',[0.05 0.1 0.85 0.8]);
        pcolorjw(nlon,nlat,nmcsst);
        caxis(clrlim);
        
        hold on;
        [cc,hh] = contour(blon,blat,bdpth,[40:10:100 200 500 1000 2000 3000]);
        set(hh,'color','k','linewidth',1);
        [cc,hh] = contour(blon,blat,bdpth,[100 200]);
        set(hh,'color','k','linewidth',2);
        clabel(cc,hh,'color',[0.5 0.5 0.5],'fontsize',15,'labelspacing',350);
        pnc;
        pmc('fill',[0.5 0.5 0.5],[0.5 0.5 0.5]);
        xlim(plonlim);ylim(platlim);
        set(gca,'xtick',[plonlim(1):0.5:plonlim(2)],'ytick',[platlim(1):0.25:platlim(2)]);
        set(gca,'layer','top','tickdir','out','fontsize',15);
        %ylabel('Latitude');
        hclb = colorbar;
        set(hclb,'location','eastoutside');
        set(hclb,'pos',[0.75 0.2 0.02 0.6],'ytick',[clrlim(1):2:clrlim(2)],'tickdir','out','fontsize',15);
        text(plonlim(2)+0.1,41,'(^oC)','fontsize',15);
        %xlabel('Longitude');
        plot(mrlon,mrlat,'pk','markersize',6,'linewidth',2);
        plot(zdlon,zdlat,'dk','markersize',6,'linewidth',2);
                
        if addadcp
            adcp   = load(adcpfile);
            trind  = find(adcp.adcptime > tstime(tranum(1)) & adcp.adcptime < tstime(tranum(end)+1));
            tru    = adcp.u(:,trind);
            trv    = adcp.v(:,trind);
            trlon  = adcp.adcplon(trind);
            trlat  = adcp.adcplat(trind);
            trtime = adcp.adcptime(trind);
            surfu  = tru(1,:);
            surfv  = trv(1,:);
            hqv    = quiver(trlon(:),trlat(:),surfu(:).*quvscale,surfv(:).*quvscale);
            set(hqv,'autoscale','off','color','k','linewidth',2);
            hqv2   = quiver(trlon(:),trlat(:),surfu(:).*quvscale,surfv(:).*quvscale);
            set(hqv2,'autoscale','off','color','w','linewidth',0.5);
            hqvs   = quiver(-70.1,39.4,uvscale.*quvscale,0);
            set(hqvs,'autoscale','off','color','k','linewidth',2);
            hqvs2   = quiver(-70.1,39.4,uvscale.*quvscale,0);
            set(hqv2,'autoscale','off','color','w','linewidth',0.5);
            text(-70.1,39.3,[num2str(uvscale) ' m s^-^1'],'fontsize',15,'color','k');
            figname = [figdir 'sst_PA_' datestr(time,'yyyymmdd') '_' datestr(time,'HHMMSS') '_adcptrsct' num2str(tranum(1)) '-' num2str(tranum(end)) '.png'];
            tttxt = ['SST on ' datestr(time) ' + ADCP at transect ' num2str(tranum(1)) '-' num2str(tranum(end))];
        else
            figname = [figdir 'sst_PA_' datestr(time,'yyyymmdd') '_' datestr(time,'HHMMSS') '.png'];
            tttxt = ['SST on ' datestr(time)];
        end
        htt = title(tttxt);
        set(htt,'fontsize',14,'fontname','Times');
        
        print('-dpng','-r300',figname);
        disp(['Printed figure file: ' figname]);
    else
        disp(['Cannot find data for ' datestr(sstime(ifile)) ' and skipping...']);
    end
    
end
