clear

%sstime(1) = datenum([2019 07 01 23 12 00]);
%sstime(1) = datenum([2019 07 04 00 23 00]);
%sstime(1) = datenum([2019 07 04 10 38 00]);
%sstime(1) = datenum([2019 07 06 00 00 00]);
sstime(1) = datenum([2019 07 08 01 10 00]);

clrlim = [15 25];

addadcp   = 1;
rmtide    = 1;  % swith to remove surface tide in adcp data
readurl   = 0;  % 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_Thompson_TN368_0703-0719/data/figs/sst/';
savemat   = 1;  % switch to save mat file
matdir    = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/sst/mat/';

if addadcp
    % ADCP data
    tranum   = [5:7];   % the transect number for adding adcp surface current
    quvscale = 0.5;   % quiver vector scale
    uvscale  = 0.2; %
    adcpfile = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/adcp/os75bb_uv_blank8mbin8m.mat';
    % read transect start time
    addpath ~/research/cruise/2019_Thompson_TN368_0703-0719/data/scripts/
    tn368_transect_startime
end

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
plonlim = [-71.85 -69.65];
platlim = [39 41.1];    % northeast

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];

% REMUS way points
auvsname = {'A18','A12','R1','R2','R3','R4','A12'};
auvlon   = [-70.830 -70.830 -70.917 -70.917 -70.745 -70.745 -70.830];
auvlat   = [ 39.620  40.010  40.010  39.880  39.880  40.010  40.010];

% streamer cross-section stations
stmrname = {'S1',    'S2',    'S3',    'S4',    'S5',    'S6',    'S7',    'S8'};
stmrlon  = [-71.0000 -70.8833 -70.7667 -70.6500 -70.5333 -70.4167 -70.3000 -70.1833];
stmrlat  = [ 39.25    39.25    39.25    39.25    39.25    39.25    39.25    39.25  ];

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
    matfiles = dir(matdir);
    for iff = 3:length(matfiles)
        matfile = [matfiles(iff).folder '/' matfiles(iff).name];
        tmptime = load(matfile,'time');
        udtime(iff-2) = tmptime.time;
    end
    tmplon  = load(matfile,'lon');
    tmplat  = load(matfile,'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
            matfile = [matfiles(tind+2).folder '/' matfiles(tind+2).name];
            tmcsst  = load(matfile,'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));
        
        if contains(matfile,'20190522_1241')
            % to manually mitigate the lat-lon offset
            nlat = nlat + 0.12;
            nlon = nlon - 0.04;
            warning(['Manually shifted lat-lon to correct the projection error!']);
        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);
        plot(auvlon,auvlat,'+k','markersize',6,'linewidth',2);
        plot(stmrlon,stmrlat,'ok','markersize',6,'linewidth',2);
                
        if addadcp
            gpsfile = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/gps/gps.mat';
            gps     = load(gpsfile);
            gpslon2 = gps.lon(1:60:end); % decimated gps for computing average speed
            gpslat2 = gps.lat(1:60:end);
            gpstim2 = datenum(gps.gpsdate(1:60:end));
            dglon2  = (gpslon2(2:end)-gpslon2(1:end-1))./180.*pi;
            dglat2  = (gpslat2(2:end)-gpslat2(1:end-1))./180.*pi;
            dgtime2 = (gpstim2(2:end)-gpstim2(1:end-1)).*86400;
            dgdist2 = sqrt(dglon2.^2 + dglat2.^2) .* rearth;
            gspeed2 = dgdist2 ./ dgtime2;
            gpstimc = (gpstim2(1:end-1) + gpstim2(2:end))./2;            
            
            adcp   = load(adcpfile);
            trind  = find(adcp.adcptime > tstime(tranum(1)) & adcp.adcptime < tstime(tranum(end)+1));
            tru    = adcp.u(:,trind);
            trv    = adcp.v(:,trind);
            trlon0 = adcp.adcplon(trind);
            trlat0 = adcp.adcplat(trind);
            trtime = adcp.adcptime(trind);
            % surfu0 = tru(1,:);
            % surfv0 = trv(1,:);
            surfu0 = tru(2,:);
            surfv0 = trv(2,:);
                        
            trsped = interp1(gpstimc,gspeed2,trtime);
            smvind = find(trsped > 1);
            surfu  = surfu0(smvind);
            surfv  = surfv0(smvind);
            trlon  = trlon0(smvind);
            trlat  = trlat0(smvind);
            trtime1= trtime(smvind);
            if rmtide
                [utide,vtide] = otps_tidecurrent_function(trlon,trlat,trtime1);
                surfu = surfu - utide;
                surfv = surfv - vtide;
            end
            
            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 'satOOI_sst_' 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))];
            if rmtide
                figname = [figname(1:end-4) '_detided.png'];
                tttxt = [tttxt ' (detided)'];
            end
        else
            figname = [figdir 'satOOI_sst_' 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
