% data_available.m  An mfile to generate a plot of data availability for
%                   the South Carolina Coastal Erosion Study data report.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
% Sullivan, C.M., Warner, J.C., Martini, M.A., Voulgaris, G., 
% Work, P.A., Haas, K.A., and Hanes, D.H. (2006) 
% South Carolina Coastal Erosion Study Data Report for Observations
% October 2003 - April 2004., USGS Open-File Report 2005-1429.
%
% Program written in Matlab v7.1.0 SP3
% Program ran on PC with Windows XP Professional OS.
%
% "Although this program has been used by the USGS, no warranty, 
% expressed or implied, is made by the USGS or the United States 
% Government as to the accuracy and functioning of the program 
% and related program material nor shall the fact of distribution 
% constitute any such warranty, and no responsibility is assumed 
% by the USGS in connection therewith."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

% C. Sullivan,  05/30/06,   version 1.0
% This code gets data availability information from all the BBV data files.
% There are 7 panels: 1) water flow obs., 2) surface wave obs. 3) adv/pcadp
% water flow obs., 4) pressure obs., 5) CT obs., 6) sediment obs., and 7)
% sonar obs.  

% Set the color scheme
siteClr = nan(9,3);
col=get(0,'defaultaxescolororder');
siteClr(1:4,:) = col(1:4,:);
siteClr(5,:) = [.9 .6 .2]  ;
siteClr(6:7,:) = col(5:6,:);
siteClr(8,:)=[.4 .4 .4];
siteClr(9,:)=[0 0 0];

% Set x-limits and ticks.  X-limits extend 5% beyond startTime and stopTime
startTime=julian(2003,10,25,00);                   
stopTime=julian(2004,04,25,00);                    
timediff=stopTime-startTime;                               
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; 
xt=[];
for yy=2003:2004
    for mm=1:12
        tic=julian(yy,mm,01,00); 
        if tic>=xlims(1) & tic<=xlims(2)
            xt=[xt; tic]; 
        end
    end
end

% Default time arrays for first and second deployment.  Times are based on
% deployment and recovery dates from the mooring log
d1Time = julian(2003,10,28):15/(60*24):julian(2004,1,23);
d2Time = julian(2004,1,29):15/(60*24):julian(2004,4,24);
d1Data = ones(size(d1Time));
d2Data = ones(size(d2Time));

% Create the figure and axes
figure
dolandscape
fixpaper
adcpAx=subplot(7,1,1,'align');
hold on
waveAx=subplot(7,1,2,'align');
hold on
advAx=subplot(7,1,3,'align');
hold on
presAx=subplot(7,1,4,'align');
hold on
salAx=subplot(7,1,5,'align');
hold on
sedAx=subplot(7,1,6,'align');
hold on
sonarAx=subplot(7,1,7,'align');
hold on

%%%%%%%%%%%%%%%%%%%%%% Water Flow Observations %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% from adcps, the argonaut and aquadopp. data from the adcp at site 5 could
% not be recovered due to excessive tilt of the instrument. site 7A did not
% have an adcp, but instead had advs and a pcadp. in the fall, data at site
% 4 was only recorded until mid-nov, when the data card filled up because
% of recording waves continuously.
if 1
% Define BBV netCDF file names for water flow data
adcpFiles = {'7201adc-a.nc','7401adc-a.nc',...  % site 1
             '7221adc-a.nc','7421adc-a.nc',...  % site 2
             '7241adc-a.nc','7441adc-a.nc',...  % site 3
             '7331adc-a.nc','7531adc-a.nc',...  % site 4
             ''            ,''            ,...  % site 5
             '7263arc-a.nc','7463arc-a.nc',...  % site 6
             ''            ,''            ,...  % site 7A
             '7321adc-a.nc','7521wh-a.nc',...   % site 7B
             '7311aqc-a.nc','7511aqc-a.nc'};    % site 8
nAdcpFiles = length(adcpFiles);

proxyVar = 'u_1205';
for f = 1:nAdcpFiles
    if ~isempty(adcpFiles{f});
        nc = netcdf(adcpFiles{f});
        mooring = nc.MOORING(:);
        if ~ischar(mooring), mooring = num2str(mooring);, end
        site = nc.DESCRIPT(6);
        if strcmp(site,'7'), site = '7B';, end
        if isequal(f,8), site = '4';, end
        time = nc{'time'}(:) + nc{'time2'}(:)/3600/1000/24;
        data = nc{proxyVar}(:,1); %use bottom bin
        theFillValue = nc{proxyVar}.FillValue_(:);
        badInd = find(data >= theFillValue);
        data(badInd) = nan;
        goodInd = find(~isnan(data));
        barStyle = '-';
    else
        if isequal(f,9) || isequal(f,13)
            time = d1Time;
            data = d1Data;
        elseif isequal(f,10) || isequal(f,14)
            time = d2Time;
            data = d2Data;
        end
        if ismember(f,[9 10]), site = '5'; barStyle = '-'; end
        if ismember(f,[13 14]), site = '7A'; barStyle = ':'; end
    end
    if mod(f,2)
        data(:) = f+1;
        fillTime = d1Time;
        fillData = d1Data .* f+1;
        cind = (f+1)/2;
        baseval = f;
        xtxt = julian(2003,11,1);
        ytxt = f+0.5;
    else
        data(:) = f;
        fillTime = d2Time;
        fillData = d2Data .* f;
        cind = f/2;
        baseval = f-1;
        xtxt = julian(2004,2,1);
        ytxt = f-0.5;
    end
    axes(adcpAx)
    BF(f)=patch([fillTime(1) fillTime(1) fillTime(end) fillTime(end)],...
                [baseval baseval+1 baseval+1 baseval], 'b');
    set(BF(f),'edgecolor', siteClr(cind,:), 'facecolor', 'none','linestyle',barStyle)
    if ~ismember(f,[9 10 13 14])
        B(f)=bar(time(goodInd), data(goodInd), 'basevalue', baseval);
        set(B(f),'edgecolor', siteClr(cind,:), 'facecolor', siteClr(cind,:),'linestyle',barStyle)
    end
    if f == nAdcpFiles
        set(B(f),'showbaseline','off')
        yl=ylabel({'Water','Flow','Observations'});
        set(yl,'position',[julian(2003,10,8) 8 0])
        set(gca,'ylim',[0 19],'ytick',[],'yticklabel',[])
        set(gca,'xlim',xlims,'xtick',xt,'xticklabel',[],'TickLength',[.003 .025])
        set(gca,'box','on','layer','top')
        [LEG,OBJ]=legend(BF([2:2:nAdcpFiles]),'Site 1','Site 2','Site 3',...
                                              'Site 4','Site 5','Site 6',...
                                              'Site 7A','Site 7B','Site 8');
        set(LEG,'Location','NorthOutside')
        set(LEG,'Orientation','Horizontal')
        figCol=get(gcf,'color');
        set(LEG,'color',figCol)
        set(LEG,'edgecolor',figCol)
        pos=get(LEG,'position');
        pobj=findobj(OBJ,'type','patch');
        for p=1:length(pobj)
            set(pobj(p),'facecolor',siteClr(p,:))
        end
        set(gcf,'renderer','painters')
        suptitle({'South Carolina Coastal Erosion Study','Data Coverage'})
        set(LEG,'position',[pos(1) .87 pos(3) pos(4)])
    end
end
 
clear f B BF LEG pobj p

end

%%%%%%%%%%%%%%%%%%%%%% Wave Observations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% waves weren't reorded at site 3 due to an instrument programming error.
% adcp at site 5 was too significantly tilted to recover data. waves
% weren't recorded at site 7A where there was only advs and a pcadp.
if 1
waveFiles={'7201adwp-cal.nc','7401adwp-cal.nc',...  % site 1
           '7221adwp-cal.nc','7421adwp-cal.nc',...  % site 2
           ''               ,'7441adwp-cal.nc',...  % site 3
           '7331adwp-cal.nc','7531adwp-cal.nc',...  % site 4
           ''               ,''               ,...  % site 5
           '7263arwp-cal.nc','7463arwp-cal.nc',...  % site 6
           ''               ,''               ,...  % site 7A
           '7321adwp-cal.nc','7521adwp-cal.nc',...  % site 7B
           '7311aqwp-cal.nc','7511aqwp-cal.nc'};    % site 8
nWaveFiles=length(waveFiles);

proxyVar = 'wh_4061';
for f = 1:nWaveFiles
    if ~isempty(waveFiles{f});
        nc = netcdf(waveFiles{f});
        time = nc{'time'}(:) + nc{'time2'}(:)/3600/1000/24;
        data = nc{proxyVar}(:); %use bottom bin
        theFillValue = nc{proxyVar}.FillValue_(:);
        %out of water data not truncated by wavesmon (or waves toolbox!)
        if mod(f,2)
            tidx1 = find(time>=julian(2003,10,28),1,'first');
            tidx2 = find(time<=julian(2004,1,23),1,'last');
        else
            tidx1 = find(time>=julian(2004,1,29),1,'first');
            tidx2 = find(time<=julian(2004,4,24),1,'last');
        end
        time = time(tidx1:tidx2);
        data = data(tidx1:tidx2);
        badInd = find(data >= theFillValue);
        data(badInd) = nan;
        goodInd = find(~isnan(data));
        barStyle = '-';
    else
        if isequal(f,5) || isequal(f,9) || isequal(f,13)
            time = d1Time;
            data = d1Data;
        elseif isequal(f,10) || isequal(f,14)
            time = d2Time;
            data = d2Data;
        end
        if ismember(f,[9 10]), barStyle = '-'; end
        if ismember(f,[13 14]), barStyle = ':'; end
    end
    if mod(f,2)
        data(:) = f+1;
        fillTime = d1Time;
        fillData = d1Data .* f+1;
        cind = (f+1)/2;
        baseval = f;
        xtxt = julian(2003,11,1);
        ytxt = f+0.5;
    else
        data(:) = f;
        fillTime = d2Time;
        fillData = d2Data .* f;
        cind = f/2;
        baseval = f-1;
        xtxt = julian(2004,2,1);
        ytxt = f-0.5;
    end
    axes(waveAx)
    BF(f)=patch([fillTime(1) fillTime(1) fillTime(end) fillTime(end)],...
                [baseval baseval+1 baseval+1 baseval], 'b');
    set(BF(f),'edgecolor', siteClr(cind,:), 'facecolor', 'none','linestyle',barStyle)
    if ~ismember(f,[5 9 10 13 14])
        B(f)=bar(time(goodInd), data(goodInd), 'basevalue', baseval);
        set(B(f),'edgecolor', siteClr(cind,:), 'facecolor', siteClr(cind,:),'linestyle',barStyle)
    end
    if f == nWaveFiles
        set(B(f),'showbaseline','off')
        yl=ylabel({'Surface','Wave','Observations'});
        set(yl,'position',[julian(2003,10,8) 8 0])
        set(gca,'ylim',[0 19],'ytick',[],'yticklabel',[])
        set(gca,'xlim',xlims,'xtick',xt,'xticklabel',[],'TickLength',[.005 .025])
        set(gca,'box','on','layer','top')
        set(gcf,'renderer','painters')
    end
end

end

%%%%%%%%%%%%%%%%%%%%%% ADV/PCADP Water Flow Observations %%%%%%%%%%%%%%%%%%
% All adv and pcadp data collected during the fall deployment is bad after
% ~20 Nov 2003 due to significant biofouling.  Only 
if 1
advFiles={''                 ,''                 ,...                %site 1
    ''                 ,''                 ,...                      %site 2
    ''                 ,''                 ,...                      %site 3
    ''                 ,''                 ,...                      %site 4
    ''                 ,''                 ,...                      %site 5
    {'adv7261vp-cal2.nc','adv7262vp-cal2.nc'},...                    %site 6 dep 1
    {'adv7461vp-cal2.nc','adv7462vp-cal2.nc'},...                    %site 6 dep 2
    {'adv7281vp-cal.nc','adv7282vp-cal2.nc','728pcvp-cal5.nc'},...   %site 7A dep 1
    {'adv7481vp-cal.nc','adv7482vp-cal2.nc','748pcvp-cal6.nc'},...   %site 7A dep 2
    ''                 ,''                 ,...                      %site 7B
    ''                 ,''                  };                       %site 8
nAdvFiles=length(advFiles);

proxyVar = 'u_1205';
for f = 1:nAdvFiles
    if ~isempty(advFiles{f});
        if length(advFiles{f}) > 1
            for n = 1:length(advFiles{f})
                nc = netcdf(advFiles{f}{n});
                time = nc{'time'}(:) + nc{'time2'}(:)/3600/1000/24;
                data = nc{proxyVar}(:,1); %use bottom bin
                theFillValue = nc{proxyVar}.FillValue_(:);
                badInd = find(data >= theFillValue);
                data(badInd)=nan;
                goodInd = find(~isnan(data));
                barStyle = '-';
                if mod(f,2)
                    if n == 1
                        data(goodInd) = f+1-.75;
                        fillTime = d1Time;
                        fillData = d1Data .* f+1-.75;
                        cind = (f+1)/2;
                        baseval = f;
                        xtxt = julian(2003,11,1);
                        ytxt = f;
                        mytxt= 'ADV1';
                    elseif n == 2
                        data(goodInd) = f+1;
                        baseval = f+.75;
                        fillData = d1Data .* f+1;
                        ytxt = f+.75;
                        mytxt= 'ADV2';
                    elseif f == 13 && n == 3 
                        data(goodInd) = f+1+.75;
                        baseval = f+1+.5;
                        fillData = d1Data .* f+1+.25;
                        ytxt = f+1+.5;
                        mytxt= 'PCADP';
                    end
                else
                    if n == 1
                        data(goodInd) = f-.75;
                        fillTime = d2Time;
                        fillData = d2Data .* f-.75;
                        cind = f/2;
                        baseval = f-1;
                    elseif n == 2
                        data(goodInd) = f;
                        baseval = f-.25;
                        fillData = d2Data .* f;
                     elseif f == 14 && n == 3
                        data(goodInd) = f+.75;
                        baseval = f+.5;
                        fillData = d2Data .* f+.75;
                    end
                end
                
                axes(advAx)
                BF(f)=patch([fillTime(1) fillTime(1) fillTime(end) fillTime(end)],...
                    [baseval baseval+.25 baseval+.25 baseval], 'b');
                set(BF(f),'edgecolor', siteClr(cind,:), 'facecolor', 'none','linestyle',barStyle)
                if ismember(f,[11:14])
                    B(f)=bar(time(goodInd), data(goodInd), 'basevalue', baseval);
                    set(B(f),'edgecolor', siteClr(cind,:), 'facecolor', siteClr(cind,:),'linestyle',barStyle)
                    set(B(f),'showbaseline','off')
                    text(xtxt,ytxt,mytxt,'fontsize',8)
                end
            end
        end
    end
    if f == nAdvFiles
        yl=ylabel({'ADV/PCADP','Water Flow','Observations'});
        set(yl,'position',[julian(2003,10,8) 12.5 0])
        set(gca,'ylim',[10.75 15],'ytick',[],'yticklabel',[])
        set(gca,'xlim',xlims,'xtick',xt,'xticklabel',[],'TickLength',[.003 .025])
        set(gca,'box','on','layer','top')
    end
end


end


%%%%%%%%%%%%%%%%%%%%%% Pressure Observations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% use adcp pressure for all except at site 7A where there were adv's
% instead of an adcp. in the fall, the pressure xducer at site 1 was
% broken, and the presssure from adcps at sites 3 and 7B was bad. site 3
% however, did have another external pressure sensor with good data.
if 1
pressureFiles={ '7201adc-a.nc','7401adc-a.nc',...        %site 1
                '7221adc-a.nc','7421adc-a.nc',...        %site 2
                '7241adc-a.nc','7441adc-a.nc'...         %site 3
                '7331adc-a.nc','7531adc-a.nc',...        %site 4
                ''            ,''           ,...         %site 5
                '7263arc-a.nc','7463arc-a.nc',...        %site 6
                'adv7281vp-cal.nc','adv7481vp-cal.nc',...%site 7A
                '7321adc-a.nc','7521wh-a.nc',...         %site 7B
                '7311aqc-a.nc','7511aqc-a.nc' };         %site 8
nPressureFiles=length(pressureFiles);

for f = 1:nPressureFiles
    if ~isempty(pressureFiles{f});
        proxyVar = 'P_4';
        if f == 13 || f == 14
            proxyVar = 'P_4023';
        end
        nc = netcdf(pressureFiles{f});
        time = nc{'time'}(:) + nc{'time2'}(:)/3600/1000/24;
        data = nc{proxyVar}(:); 
        theFillValue = nc{proxyVar}.FillValue_(:);
        badInd = find(data >= theFillValue);
        data(badInd) = nan;
        goodInd = find(~isnan(data));
        barStyle = '-';
    else
        if f == 9
            time = d1Time;
            data = d1Data;
        elseif f == 10
            time = d2Time;
            data = d2Data;
        end
        barStyle = ':';
        if f == 1 || f == 15
            barStyle = '-';
        end
    end
    if mod(f,2)
        data(goodInd) = f+1;
        fillTime = d1Time;
        fillData = d1Data .* f+1;
        cind = (f+1)/2;
        baseval = f;
        xtxt = julian(2003,11,1);
        ytxt = f+0.5;
    else
        data(goodInd) = f;
        fillTime = d2Time;
        fillData = d2Data .* f;
        cind = f/2;
        baseval = f-1;
        xtxt = julian(2004,2,1);
        ytxt = f-0.5;
    end
    axes(presAx)
    BF(f)=patch([fillTime(1) fillTime(1) fillTime(end) fillTime(end)],...
                [baseval baseval+1 baseval+1 baseval], 'b');
    set(BF(f),'edgecolor', siteClr(cind,:), 'facecolor', 'none','linestyle',barStyle)
    if ~ismember(f,[1 5 9 10 15])
        set(BF(f),'linestyle','-')
        B(f)=bar(time(goodInd), data(goodInd), 'basevalue', baseval);
        set(B(f),'edgecolor', siteClr(cind,:), 'facecolor', siteClr(cind,:),'linestyle',barStyle)
        set(B(f),'showbaseline','off')
    end
    if f == nPressureFiles
        yl=ylabel({'Pressure','Observations'});
        set(yl,'position',[julian(2003,10,8) 8 0])
        set(gca,'ylim',[0 19],'ytick',[],'yticklabel',[])
        set(gca,'xlim',xlims,'xtick',xt,'xticklabel',[],'TickLength',[.005 .025])
        set(gca,'box','on','layer','top')
        set(gcf,'renderer','painters')
    end
end

end

%%%%%%%%%%%%%%%%%%%%%% CT Observations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only show bars for sites that had CTs.  Do not include bars for all sites
% in order to retain panel clarity.
if 1
ctFiles={''           ,''           ,...                    %site 1
         {'7222sc-a.nc','7231mc-a.nc'},...                  %site 2 bottom, surface spring
         {'7422sc-a.nc','7431mc-a.nc'},...                  %site 2 bottom, surface fall
         {'7242sc-a.nc','7251mc-a.nc'},...                  %site 3 bottom, surface spring
         {'7442sc-a.nc','7451mc-a.nc'},...                  %site 3 bottom, surface fall
         ''           ,''           ,...                    %site 4
         ''           ,''           ,...                    %site 5
         {'7265sc-a.nc'},{'7465Asc-a.nc';'7465Bsc-a.nc'},...%site 6 bottom
         {'7291mc-a.nc'},{'7491mc-a.nc'},...                %site 7A surface
         ''           ,''           ,...                    %site 7B
         ''           ,''           };                      %site 8
nCtFiles=length(ctFiles);

proxyVar = 'S_40';
for f = 1:nCtFiles
    if ~isempty(ctFiles{f});
            for n = 1:length(ctFiles{f})
                if f == 12
                    nc1 = netcdf(ctFiles{f}{1});
                    nc2 = netcdf(ctFiles{f}{2});
                    time = [nc1{'time'}(:) + nc1{'time2'}(:)/3600/1000/24;...
                            nc2{'time'}(:) + nc2{'time2'}(:)/3600/1000/24];
                    data = [nc1{proxyVar}(:);nc2{proxyVar}(:)];
                    theFillValue = nc1{proxyVar}.FillValue_(:);                  
                else
                    nc = netcdf(ctFiles{f}{n});
                    time = nc{'time'}(:) + nc{'time2'}(:)/3600/1000/24;
                    data = nc{proxyVar}(:); %use middle bin
                    theFillValue = nc{proxyVar}.FillValue_(:);
                end
                badInd = find(data >= theFillValue);
                data(badInd) = nan;
                goodInd = find(~isnan(data));
                barStyle = '-';
                if mod(f,2)
                    if n == 1
                        data(goodInd) = f+.25;
                        fillTime = d1Time;
                        fillData = d1Data .* f+.25;
                        cind = (f+1)/2;
                        baseval = f;
                        xtxt = julian(2003,10,18);
                        ytxt = f+.125;
                        mytxt= 'bottom';
                        if f>=11
                            data(goodInd) = f-4+.25;
                            fillData = d1Data .* f-4+.25;
                            baseval = f-4;
                            ytxt = f-4;
                        end
                    elseif n == 2
                        data(goodInd) = f+1;
                        baseval = f+.75;
                        fillData = d1Data .* f+1;
                        ytxt = f+.75+.25;
                        mytxt= 'surface';
                        if f>=11
                            data(goodInd) = f+1-4;
                            fillData = d1Data .* f+1-4;
                            baseval = f+.75-4;
                            ytxt = ytxt-4;
                        end
                    end
                else
                    if n == 1 || f==12
                        data(goodInd) = f-.75;
                        fillTime = d2Time;
                        fillData = d2Data .* f-.75;
                        cind = f/2;
                        baseval = f-1;
                        if f>=11
                            data(goodInd) = data(goodInd)-4;
                            fillData = fillData-4;
                            baseval = baseval-4;
                        end
                    elseif n == 2
                        data(goodInd) = f;
                        baseval = f-.25;
                        fillData = d2Data .* f;
                        if f>=11
                            data(goodInd) = data(goodInd)-4;
                            fillData = fillData-4;
                            baseval = baseval-4;
                        end
                    end
                end
                
                axes(salAx)
                BF(f)=patch([fillTime(1) fillTime(1) fillTime(end) fillTime(end)],...
                    [baseval baseval+.25 baseval+.25 baseval], 'b');
                set(BF(f),'edgecolor', siteClr(cind,:), 'facecolor', 'none','linestyle',barStyle)
                if ismember(f,[3:6 11:14])
                    B(f)=bar(time(goodInd), data(goodInd), 'basevalue', baseval);
                    set(B(f),'edgecolor', siteClr(cind,:), 'facecolor', siteClr(cind,:),'linestyle',barStyle)
                    set(B(f),'showbaseline','off')
                    if f == 14
                        obj=findobj('type','text','position',[xtxt ytxt 0]);
                        delete(obj)
                        text(xtxt,ytxt,'surface','fontsize',8,'verticalalignment','middle')
                    else
                        text(xtxt,ytxt,mytxt,'fontsize',8,'verticalalignment','middle')
                    end
                end
            end
    end
    if f == nCtFiles
        yl=ylabel({'Temperature','and Salinity ','Observations'});
        set(yl,'position',[julian(2003,10,8) 6 0])
        set(gca,'ylim',[2.5 9.75],'ytick',[],'yticklabel',[])
        set(gca,'xlim',xlims,'xtick',xt,'xticklabel',[],'TickLength',[.003 .025])
        set(gca,'box','on','layer','top')
    end
end

end



%%%%%%%%%%%%%%%%%%%%%% ABS Sediment Concentration Observations %%%%%%%%%%%%
% Use ABS data only.  OBS data from both deployments was bad.
if 1
absFiles={''               ,''               ,... %site 1
          ''               ,''               ,... %site 2
          ''               ,''               ,... %site 3
          ''               ,''               ,... %site 4
          ''               ,''               ,... %site 5
          '7264abss-cal.nc','7464abss-cal.nc',... %site 6
          ''               ,'7484abss-cal.nc',... %site 7A
          '7322abss-cal.nc','7522abss-cal.nc',... %site 7B
          ''               ,''                };  %site 8
nAbsFiles=length(absFiles);

proxyVar = 'abs_trans1_mean';
for f = 1:nAbsFiles
    if ~isempty(absFiles{f});
        nc = netcdf(absFiles{f});
        time = nc{'time'}(:) + nc{'time2'}(:)/3600/1000/24;
        distFromSensor = nc{'r'}(:,1);
        nBins = length(distFromSensor);
        data = nc{proxyVar}(:,nBins/2); %use middle bin
        theFillValue = nc{proxyVar}.FillValue_(:);
        badInd = find(data >= theFillValue);
        data(badInd) = nan;
        goodInd = find(~isnan(data));
        barStyle = '-';
    else
        if ismember(f,[1 3 5 7 9 13 17])
            time = d1Time;
            data = d1Data;
        elseif ismember(f,[2 4 6 8 10 18])
            time = d2Time;
            data = d2Data;
        end
        barStyle = ':';
    end
    if mod(f,2)
        data(:) = f+1;
        fillTime = d1Time;
        fillData = d1Data .* f+1;
        cind = (f+1)/2;
        baseval = f;
        xtxt = julian(2003,11,1);
        ytxt = f+0.5;
    else
        data(:) = f;
        fillTime = d2Time;
        fillData = d2Data .* f;
        cind = f/2;
        baseval = f-1;
        xtxt = julian(2004,2,1);
        ytxt = f-0.5;
    end
    axes(sedAx)
    BF(f)=patch([fillTime(1) fillTime(1) fillTime(end) fillTime(end)],...
                [baseval baseval+1 baseval+1 baseval], 'b');
    set(BF(f),'edgecolor', siteClr(cind,:), 'facecolor', 'none','linestyle',barStyle)
    if ismember(f,[11 12 14 15 16])
        set(BF(f),'linestyle','-')
        B(f)=bar(time(goodInd), data(goodInd), 'basevalue', baseval);
        set(B(f),'edgecolor', siteClr(cind,:), 'facecolor', siteClr(cind,:),'linestyle',barStyle)
        set(B(f),'showbaseline','off')
    end
    if f == nAbsFiles
        yl=ylabel({'Sediment','Observations'});
        set(yl,'position',[julian(2003,10,8) 8 0])
        set(gca,'ylim',[0 19],'ytick',[],'yticklabel',[])
        set(gca,'xlim',xlims,'xtick',xt,'xticklabel',[],'TickLength',[.005 .025])
        set(gca,'box','on','layer','top')
        set(gcf,'renderer','painters')
    end
end

end
    
    
%%%%%%%%%%%%%%%%%%%%%% Sonar Observations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use the dates of the first and last images from the animations to
% determine where there is/isn't data. Sonar data only at site 7B.  The
% sonar at site 6 was tilted and/or flooded.
if 1
sonarFiles={''                                 ,''                                 ,... % site 1
            ''                                 ,''                                 ,... % site 2
            ''                                 ,''                                 ,... % site 3
            ''                                 ,''                                 ,... % site 4
            ''                                 ,''                                 ,... % site 5
            ''                                 ,''                                 ,... % site 6
            ''                                 ,''                                 ,... % site 7A
            {'7323Aim-cal.nc','7323Bim-cal.nc'},{'7523Aim-cal.nc','7523Bim-cal.nc'},... % site 7B
             ''                                ,''                                  };  % site 8
nSonarFiles = length(sonarFiles);

d1FirstImage=julian([2003,10,28,14,20,00]);
d1LastImage=julian([2004,1,21,20,50,00]);
d2FirstImage=julian([2004,1,29,19,50,00]);
d2LastImage=julian([2004,3,18,19,50,00]);
site = '7B'
proxyVar = 'sonar_image';
for f = 1:nSonarFiles
    if f == 15
        time = d1FirstImage:d1LastImage;
        data = ones(size(time)) .* f+1;
    elseif f == 16
        time = d2FirstImage:d2LastImage;
        data = ones(size(time)) .* f;
    end
    if mod(f,2)
        fillTime = d1Time;
        fillData = d1Data .* f+1;
        cind = (f+1)/2;
        baseval = f;
        xtxt = julian(2003,11,1);
        ytxt = f+0.5;
    else
        fillTime = d2Time;
        fillData = d2Data .* f;
        cind = f/2;
        baseval = f-1;
        xtxt = julian(2004,2,1);
        ytxt = f-0.5;
    end
    axes(sonarAx)
    BF(f)=patch([fillTime(1) fillTime(1) fillTime(end) fillTime(end)],...
                [baseval baseval+1 baseval+1 baseval], 'b');
    set(BF(f),'edgecolor', siteClr(cind,:), 'facecolor', 'none','linestyle',':')
    if ismember(f,[15 16])
        set(BF(f),'linestyle','-')
        B(f)=bar(time, data, 'basevalue', baseval);
        set(B(f),'edgecolor', siteClr(cind,:), 'facecolor', siteClr(cind,:),'linestyle','-')
        set(B(f),'showbaseline','off')
    end
    if f == nSonarFiles
        yl=ylabel({'Sea Floor','Bedform','Observations'});
        set(yl,'position',[julian(2003,10,8) 8 0])
        set(gca,'ylim',[0 19],'ytick',[],'yticklabel',[])
        set(gca,'xlim',xlims,'xtick',xt,'TickLength',[.005 .025],'xticklabel',[])
        set(gca,'box','on','layer','top')
        for tic=1:length(xt)
            xtg=gregorian(xt(tic));
            xtd=datenum(xtg);
            xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)};
            text(xt(tic),-1,xtl,'horizontalalignment','center',...
                'verticalalignment','top')
        end
        set(gcf,'renderer','painters')
    end
end

end