Managing and Plotting Buoy Data in MATLAB

In this page, the approach taken for building a plot of an average yearly wave climate for the Umpqua Offshore buoy is described. The two parameters pulled from the publicly available data sets are significant wave height, Hs, and dominant wave period, Td. A description of these variables can be found on the NDBC site here.

The purpose of this page is to provide m code that creates plots to disseminate a lot of information clearly and concisly. For example, we will be using yearly parameter data from years 2005 - 2010 (download .txt files here) to can create a two-dimensional plot that looks something like Fig. 1.

example figure

Figure 1: An example plot of a wave climate without proper formatting.


To begin, first download MATLAB's hdrload.m as we will be using it to parse the text files. Next, implement the following code, where the input is the raw text files downloaded from the NDBC website and the output is an x by y matrix for each specified parameter (x [col] = yearly data and y [row] = number of parameter entries).

%This script opens the raw .txt files downloaded from NDBC and writes the
%wvht, dpd, and date to specific variables. Where each column represents a
%year (col 1 = 2005, ... col 6 = 2010) and so on. 

%% - Load the .txt files from the NDBC website
%First, create a cell matrix that will hold all of the data, from here we
%can pull specific parameters (wave height, wave period, etc.)
data = cell(6,1);

%Next, use the MATLAB function "hdrload" to extract the header and data
[h data{1,1}] = hdrload('umpqua2005.txt');
[h data{2}] = hdrload('umpqua2006.txt');
[h data{3}] = hdrload('umpqua2007.txt');
[h data{4}] = hdrload('umpqua2008.txt');
[h data{5}] = hdrload('umpqua2009.txt');
[h data{6}] = hdrload('umpqua2010.txt');

%In this case, there are six years of data to aggregate
for n = 1:6
    
    %Create and assign data to intermediate variables. The entire column of
    %data (9 = Hs, 10 = Td) is assigned to wv and dp respectively.
    wv = data{n,1}(:,9); dp = data{n,1}(:,10);
    
    %Assign wv and dp to master list of wvht and dpd. Here the intermediate
    %variables are assigned to a master matrix of parameters. Where the
    %first year of data (2005) is the first column in wvht & dpd.
    wvht(1:length(wv),n) = wv;
    dpd(1:length(dp),n) = dp;
    
    %Knowing that there are anomalies in the raw data files that were not 
    %extracted, parse every observation of the current parameter stream
    %using the variable m. 
    for m = 1:length(wv)
        
        %Assign a date to every row that contains a wvht and dpd
        da = datenum(data{n,1}(m,1), data{n,1}(m,2), data{n,1}(m,3),...
		data{n,1}(m,4), data{n,1}(m,5), 0);
        
        date(m,n) = da;
        
        %Find any obs error and replace with averages of nearest neighbors
        if wvht(m,n) == 99 %99 is the value recorded for false readings
            wvht(m,n) = (wvht(m+1,n)+wvht(m-1,n))/2;
            dpd(m,n)  = (dpd(m+1,n)+dpd(m-1,n))/2;
            
            %If there are more than one bad recording next to each other,
            %using neighbors that are farther away.
            
            if wvht(m,n) >= 25 %25 is an arbitrary valued
                wvht(m,n) = (wvht(m+2,n)+wvht(m-2,n))/2;
                dpd(m,n) = (dpd(m+2,n)+dpd(m-2,n))/2;
            end
        end
        
        %Clear da as an intermediate variable
        clear da
    end
    
    %Clear wv and dp because they are intermediate variables
    clear wv dp
end
%Clear data because it's huge
clear data

%Now, the master variables have been built and their failed recordings
%removed and replaced with averages of their nearest neighbors.
%Save the variables to disk and 
save('wvht.mat','wvht');
save('dpd.mat','dpd');
save('date.mat','date');

Download: extract_parameters.m

Now that we have two variables, wvht and dpd (both matrices ~8000 x 6), we can use them to build a histogram of the average wave climate for an average year. This next script explains how to accomplish this problem.


% This quick script searches the parameter dpd for the exact row that ends
% the reported data set. This is necessary because each year has a
% different data set length. The variable "ennd" is used to build two sets:
% one where the time step is 1 hour and the other is 0.5 hours. Sometime
% around the beginning of February 2008, the reports started being recorded
% every half hour.

for s = 1:6
    ennd(s) = length(find(dpd(:,s)>0));
end

%Build wave climate from records spaced 1 hour apart, note row number 697
%in year 2008 is the last hour recording
wave.a = [wvht(1:ennd(1),1); wvht(1:ennd(2),2); wvht(1:ennd(3),3); wvht(1:697,4)];
period.a = [dpd(1:ennd(1),1); dpd(1:ennd(2),2); dpd(1:ennd(3),3); dpd(1:697,4)];
sea_state.a = [period.a wave.a];
hours.a = sum(ennd(1:3))+697;

%Build wave climate from records spaced 0.5 hours apart
wave.b = [wvht(698:ennd(4),4); wvht(1:ennd(5),5); wvht(1:ennd(6),6)];
period.b = [dpd(698:ennd(4),4); dpd(1:ennd(5),5); dpd(1:ennd(6),6)];
sea_state.b = [period.b wave.b];
hours.b = sum(ennd(5:6))+(ennd(4)-698); 
%Set total half hours reported in this climate to whole hours
hours.b = hours.b/2;

%Net time - not exaclty 6 years
hours.c = (hours.a + hours.b)/24/365;

%% Build a three dimensional histogram using the function "hist3"

% Set the boundaries for the historgram
dx1 = 1; % wave period distributed in bins of 1 sec 
dy1 = 0.5; % wave height distributed in bins of 0.5 m

% Set the center of the bins,
x1 = 1.5:dx1:22.5; % wave period
y1 = 0.25:dy1:10.75; % wave height

% This variable will be used as the framework for the histogram
centers = {x1 y1};

% Extract the histogram data
num.a = hist3(sea_state.a, 'Ctrs',centers);
num.b = hist3(sea_state.b, 'Ctrs',centers);

% Divide num.b by 2 to set the count values to hours
num.b = num.b/2;

% Add both matrices together and divide by the total number of hours the
% data set represents. This variable is a 22 x 22 matrix
num.c = (num.a + num.b)/hours.c;

%% Plot the histogram data

% Take the inverse of the num.c matrix. This puts the wave period on the x-axis
% and the wave height on the y-axis.
plot_hours = num.c';
% Add a blank row around the outer edges to make it possible to see all the
% entries on the plot
plot_hours( size(num.c,1) + 1 ,size(num.c,2) + 1 ) = 0;


%create figure object to assign the wave climate
figure

% This function actually plots a pseudocolor plot with the standard
% colormap from a 23 x 23 matrix. 
h = pcolor(plot_hours);

% Set the colormap to gray and inverse the colors: black for the highest
% values and white for the lowest.
colormap(gray)
colormap(flipud(gray))


%% - Use colormapeditor to change colors (cd max = 550)

gray colormap

Figure 2: Original gray colormap.

However, if we are to overlay numbers on this plot in black, then we need to change how dark the highest values are plotted. This is done using the colormapeditor function from the command window (Fig. 3).

colormapeditor

Figure 3: The colormapeditor toolbox.

After executing this procedure, the figure should resemble Fig. 4:

gray

Figure 4: New colormap scheme.

 

Now, continue with executing the following code:

%% Countinue with editing figure

% Change the edge color to gray 
set(h,'linestyle',:,'edgecolor',[0.5 0.5 0.5])

% Overlay the number of hours for the yearly average for each bin 
[m,n]=size(plot_hours);

for r = 1:m
    for s = 1:n
        if plot_hours(r,s) > 0 && r ~= m && s ~= n
            text((s+dx1/2)-0.1,r+dy1/2,num2str(plot_hours(r,s),'%5.1f'), 'fontsize',8)
            %         elseif plot_hours(r,s) > 0 && s == n
            %             text(s,r,num2str(plot_hours(r,s)))
            %         elseif plot_hours(r,s) > 0 && r == m
            %             text(x1(r)-dx1/2,y1(s)+dy1/2,num2str(plot_hours(r,s)))
        else
        end
    end
end

% Assign axis labels to the figure
ylabel('Significant Wave Height [m]','fontsize',12,'interpreter','latex')
xlabel('Dominant Wave Period [s]','fontsize',12,'interpreter','latex')

% Apply axis changes (plot_hours matrix is in standard 1,2,3  x 1,2,3 index,
% however we need to assign the labels the correct sea state index

% Default xlim = [1 23], move plot left to gain more space to plot
xlim([2 23])
dplimit = get(gca,'xlim');
% Wave height values
yticklabel1 = {'0','1','2','3','4','5','6','7','8','9','10','>11'};
ytick1 = [1 3 5 7 9 11 13 15 17 19 21 23];
% Wave period values
xticklabel1 = {'2','4','6','8','10','12','14','16','18','20','22','>23'};
xtick1 = [2 4 6 8 10 12 14 16 18 20 22 23];
set(gca,'YTickLabel',yticklabel1, 'YTick',ytick1, 'XTickLabel',xticklabel1,...
    'XTick',xtick1);
xbarlimit = get(gca,'ylim');

% Prepare the figure for plotting on letter size paper
picture = fillPage(gcf,'margins',[0 0 0 4.5],'papersize',[11 8.5]);
print(gcf, '-dpdf','-r500', 'average_climate.pdf')
  

The final product from the previous code can be seen in Fig. 5. which shows the number of hours for an average year that the Umpqua Offshore buoy experiences a particular sea state.

 


Next, the dominant wave period histogram can be plotted.

 %% Wave period histogram 
figure
% Set the boundaries of the matrix
dp_edges = [1:22,inf];

% Similar to the previous organization with both wave parameters, there is
% discrepancy with the hour and half hour recording times
dp_indp.a = histc(period.a,dp_edges);
dp_indp.b = histc(period.b,dp_edges); dp_indp.b = dp_indp.b/2;

% Find the total hours
dp_indp.c = (dp_indp.a + dp_indp.b)/hours.c;

% Find the average wave period
meanT.a = mean(period.a); meanT.b = mean(period.b);
meanT.c = (hours.a * meanT.a + hours.b * meanT.b) /(hours.a + hours.b);

% Find the standard deviation 
stdT.a = std(period.a); stdT.b = std(period.b);
stdT.c = sqrt((hours.a * stdT.a^2 + hours.b * stdT.b^2)/(hours.a + hours.b) ...
    + (hours.a * hours.b) * ((meanT.a - meanT.b)^2 /(hours.a + hours.b)^2));

floor_meanT = floor(meanT.c);

% Actually plot the histogram
stairs(dp_indp.c,'k','LineWidth',3)
% Make background white
grid off; hold on; box off;

% line([meanT.c meanT.c],[0 1500],'r..','linewidth',3)
% text(10,10,['Avg = ', num2str(meanT.c), 'SD = ', num2str(stdT.c)],'%5.1f')
xlabel('Dominant Wave Period [sec]','fontsize',16,'interpreter','latex')
ylabel('Occurrence [hr]','fontsize',16,'interpreter','latex')

% Set axis label values and axis separations limits equal to the master
% plot in subplot
xlim(dplimit)

print(gcf, '-dpdf','-r500', 'average_period.pdf')
   

Download: plot_average_year.m

wave period

Figure 6: Histogram of dominant wave period from entire data set.

If you have any questions, send me an email.