Please support this Forum and Neven's Blog

Author Topic: Matlab code to plot IMB Buoy profiles  (Read 555 times)

seaicesailor

  • ASIF Upper Class
  • Posts: 1191
    • View Profile
Matlab code to plot IMB Buoy profiles
« on: August 01, 2016, 07:31:17 PM »
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
« Last Edit: August 01, 2016, 07:50:44 PM by seaicesailor »