Here is the Matlab code I used to create the gif of IMB Buoy profiles
- First copy the code below and put it into a text file called plotBuoyTemperatureProfile.m in whatever folder
- Download the csv file that is provided in the IMB buoy page and put it in the same folder:
http://imb.erdc.dren.mil/irid_data/2015F_clean.csv
- Open Matlab, go to the folder, and write:
[ data, filteredData ] = plotBuoyTemperatureProfile( '2015F_clean.csv', '07-15-2016', '08-01-2016', 1 );
Arguments are: csv file name, initial date, final date, and profile step. If profile step = 1, it will plot all profiles within the dates (six per day) which is very slow (see example attached). Profile step = 6 will plot one profile per day.
The outputs are the data and the filtered data array of the thermistor temperatures.
The function will automatically create a folder called "figs" and will save all the individual png of each plotted profile, as well as the gif. I set the frame delay at 0.3 seconds but it can be changed at line 83
Without Matlab, it may work with Octave but probably will need some modifications.
% -------------------------------------------------------------------------
% [ data, filteredData ] = plotBuoyTemperatureProfile(
% filename, initialDate, finalDate, profileStep )
% -------------------------------------------------------------------------
% s.i.s. July 2016
% - filename: mass balance buoy csv filename such as '2015F_clean.csv'
% - initialDate: date in string 'mm-dd-yy' format
% - finalDate: date in string 'mm-dd-yy' format
% - profileStep: plot only one per profileStep profiles
% For Buoy 2015F, there are six profiles per day.
% -------------------------------------------------------------------------
% The function generates one png plot per profile, also outputs original
% and filtered arrays of thermistors temperatures
% GPS, air temp, thickness, header, rows with no thermistor temperature
% data... not used
% -------------------------------------------------------------------------
function [ data, filteredData ] = ...
plotBuoyTemperatureProfile( filename, initialDate, finalDate, profileStep )
fid = fopen(filename, 'r');
dataRaw = textscan(fid,'%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s','Delimiter', ',' , 'Headerlines', 1);
fclose(fid);
% Parse data, remove rows with blank temperature profiles
j = 0;
for idx = 1 : numel(dataRaw{1})
if ~isempty(dataRaw{11}{idx}(2:end-1))
j = j + 1;
dateStr{j} = datestr(dataRaw{1}{idx}(2:end-1),'mm-dd-yy');
dateNum(j) = datenum(dataRaw{1}{idx}(2:end-1));
% data array with thermistor profiles
for k = 1:31
data ( j, k ) = str2double(dataRaw{k+10}{idx}(2:end-1));
end
end
end
sizeArray = size(data);
% Filter in the time direction: a parabolic window filter with size 12
% elements (two days) to remove fluctuations below the 1-day scale
aux = data;
for n = 1: sizeArray(2)
aux(:,n) = filterParabolicWindow(data(:,n), 18);
end
% Filter in the direction of the thermistors: a parabolic window filter
% with 4 elements (40 centimeters) to remove fluctuations below the 20
% centimeter scale.
filteredData = aux;
for n = 1: sizeArray(1);
filteredData(n,:) = filterParabolicWindow(aux(n,:), 4);
end
% Find the initial and final profiles depending on the input dates
initialProfile = find(dateNum >= datenum(initialDate,'mm-dd-yy'),1,'first');
finalProfile = find(dateNum >= datenum(finalDate,'mm-dd-yy'),1,'first');
% Check correct limits
if isempty(initialProfile), initialProfile = 1; end
if isempty(finalProfile), finalProfile = sizeArray(1); end
finalProfile=max(finalProfile, initialProfile);
% Finally plot the selected profiles within the dates provided
figFolder = 'figs';
mkdir(figFolder);
gifName = fullfile( figFolder, [ filename '_' initialDate '_' finalDate '.gif'] );
for iG = initialProfile:profileStep:finalProfile
pFig = plotBuoyprofile(1:sizeArray(2), filteredData(iG,1:sizeArray(2)), dateStr{iG}, iG, figFolder);
drawnow
frame = getframe(1);
im = frame2im(frame);
[imind,cm] = rgb2ind(im,256);
if iG == initialProfile;
imwrite(imind,cm,gifName,'gif', 'Loopcount',inf);
else
imwrite(imind,cm,gifName,'gif','WriteMode','append','DelayTime',0.3);
end
close(pFig);
end
end
function pfig = plotBuoyprofile(thermistorNumberArray, temperatureArray, dateStr, idx, figFolder)
% Create figure
pfig = figure('PaperSize',[20.98404194812 29.67743169791],...
'InvertHardcopy','off',...
'Color',[1 1 1]);
% Create axes
axes1 = axes('Parent',pfig,'YGrid','on','XGrid','on','FontSize',14);
xlim(axes1,[0 250]);
ylim(axes1,[-3 2]);
box(axes1,'on');
hold(axes1,'all');
% Create plot
plot(thermistorNumberArray*10, temperatureArray,'MarkerSize',4,'Marker','o','LineWidth',3,'Color',[0 0 1]);
% Create xlabel
xlabel('Depth(cm)','FontSize',14);
% Create ylabel
ylabel('Temperature (^oC)','FontSize',14);
% Create title
title([ 'Buoy 2015F profile - ' dateStr],'FontSize',14);
% Create multiple lines using matrix input to plot
plot1 = plot(thermistorNumberArray*10,thermistorNumberArray*0,'LineStyle','--');
plot2 = plot(thermistorNumberArray*10,-1.8*thermistorNumberArray./thermistorNumberArray,'LineStyle','--');
set(plot1(1),'LineWidth',3,'Color',[0 0 0]);
set(plot2(1),'LineWidth',2,'Color',[0.847058832645416 0.160784319043159 0]);
saveas(pfig, fullfile(figFolder, ['date_' dateStr '_' num2str(idx)]), 'png');
end
% -------------------------------------------------------------------------
% Classic window filter with parabolic weighting
% -------------------------------------------------------------------------
function filteredVariable = filterParabolicWindow( inputVariable, windowSize )
sizeInputVariable = size(inputVariable);
numberOfElements = numel(inputVariable);
inputVariable = reshape(inputVariable, 1, numberOfElements);
% Symmetric parabola equal to zero at borders
rangeWindow = 1:2*windowSize+1;
parabolicWindow = (rangeWindow-1).*(2*windowSize+1- rangeWindow);
filteredVariable = inputVariable;
for ifilter = 1:numberOfElements
% Window shrunk at sides of domain
indexVar = max(1,ifilter - windowSize):min(ifilter + windowSize, numberOfElements);
indexHat = indexVar - ifilter + windowSize + 1;
% Apply filtering over the window range
filteredVariable(ifilter) = sum(parabolicWindow(indexHat).*filteredVariable(indexVar) )/...
sum(parabolicWindow(indexHat));
end
filteredVariable = reshape(filteredVariable, sizeInputVariable);
end
