Source Code:

%-----------------------------------------------------------------------
% Name:
%   .dbase/web_programs/plotpro/matlab_scripts/difference.m
% Description:
%   Creates hs compendium like difference plot. A few lines are 
%   pre-pended to this script.
% Called by:
%   .f90/plot_utils.make_plot
%-----------------------------------------------------------------------

global interp_t interp_pm numdate count inc

% ADD PATH TO ACCESS MATLAB SCRIPTS

path(path,'/project/dbase/web_programs/plotpro/matlab_scripts');

wkg = wkg_path;

% GET STATION INFO

info1 = read_info_file([pid,'01'],wkg);
info2 = read_info_file([pid,'02'],wkg);


% specify the color of the second plotted line.
dg = [ 1 .6 .1 ];

% load data

eval(['pm_file1=load(''',wkg,'list.',pid,'01'');']);
eval(['pm_file2=load(''',wkg,'list.',pid,'02'');']);

figure('Visible','off','Position',[100,100,900,700],'PaperPosition',[0 0 7.75 7.75]);

pmlabel = ['Difference of Significant Wave Heights (1-2) (m)'];
pmcolumn = 6;

% read in date, parameter values

x=size(pm_file1);
yr=pm_file1(1:x(1),1);
mo=pm_file1(1:x(1),2);
dd=pm_file1(1:x(1),3);
hr=pm_file1(1:x(1),4);
mm=pm_file1(1:x(1),5);
pm=pm_file1(1:x(1),pmcolumn);

max1 = max(pm);
min1 = min(pm);

% convert time variables to datenumbers for plots and labels

numdate=datenum(yr,mo,dd,hr,mm,0);

start_time1 = numdate(1);
end_time1 = numdate(x(1));




% draw parameter line, including gaps where data is missing
% this algorithm uses a four-hour maximum interval
% between data runs; if runs are farther apart than
% 4 hours (.1667 days), the time series is considered broken and
% a new line is begun

% in addition, construct interpolated data set for comparisons

interp_t = [];
interp_pm = [];
inc = 0;
count = 1;
total = length(numdate);
while count < total
  if numdate(count+1) <= numdate(count)+(.1667)
    inc = inc + 1;
  else
    if inc > 0
      interpolate;
      hold on;
      plot(numdate(count-inc:count),pm(count-inc:count),'b-');
      inc = 0;
    end
  end
  count = count + 1;
end

interpolate;
hold on;
plot(numdate(count-inc:count),pm(count-inc:count),'b-');

times_1 = interp_t;
vals_1 = interp_pm;


% repeat above steps for second parameter file

x=size(pm_file2);
yr=pm_file2(1:x(1),1);
mo=pm_file2(1:x(1),2);
dd=pm_file2(1:x(1),3);
hr=pm_file2(1:x(1),4);
mm=pm_file2(1:x(1),5);
pm=pm_file2(1:x(1),pmcolumn);

max2 = max(pm);
min2 = min(pm);

numdate=datenum(yr,mo,dd,hr,mm,0);

start_time2 = numdate(1);
end_time2 = numdate(x(1));

% draw parameter line, including gaps where data is missing

interp_t = [];
interp_pm = [];
inc = 0;
count = 1;
total = length(numdate);
while count < total
  if numdate(count+1) <= numdate(count)+(.1667)
    inc = inc + 1;
  else
    if inc > 0
      interpolate;
      hold on;
      plot(numdate(count-inc:count),pm(count-inc:count),'Color',dg);
      inc = 0;
    end
  end
  count = count + 1;
end

interpolate;
hold on;
plot(numdate(count-inc:count),pm(count-inc:count),'Color',dg);

times_2 = interp_t;
vals_2 = interp_pm;

% calculate the difference wherever data exists for the same time
% in the two interpolated arrays

difference = [];
diff_times = [];
count1 = 1;
count2 = 1;
total1 = length(times_1);
total2 = length(times_2);
 
while (count1 < total1) & (count2 < total2) 
  if times_1(count1) == times_2(count2)
    difference = [difference vals_1(count1)-vals_2(count2)];
    diff_times = [diff_times times_1(count1)];
    count1 = count1 + 1;
    count2 = count2 + 1;
  elseif times_1(count1) > times_2(count2)
    count2 = count2 + 1;
  else
    count1 = count1 + 1;
  end
end

% plot for the ratio data, skipping over 4-hour gaps

inc = 0;
count = 1;
total = length(diff_times);
while count < total
  if diff_times(count+1) <= diff_times(count)+(.1667)
    inc = inc + 1;
  else
    if inc > 0
      gap_times = [diff_times(count-inc:count-inc) diff_times(count-inc:count) diff_times(count:count)];
      gap_difference = [0 difference(count-inc:count) 0];
      hold on;
      plot(gap_times,gap_difference,'r','markersize',4);
      inc = 0;
    end
  end
  count = count + 1;
end

gap_times = [diff_times(count-inc:count-inc) diff_times(count-inc:count) diff_times(count:count)];
gap_difference = [0 difference(count-inc:count) 0];
hold on;
plot(gap_times,gap_difference,'r','markersize',4);

start_time = min([start_time1 start_time2]);
end_time = max([end_time1 end_time2]);

timespan=end_time-start_time;
timediv=(timespan/8);
xar=[start_time:timediv:end_time];
xlab1=datestr(xar,'mm/dd/yy');
xlab2=datestr(xar,'HH:MM');

% make a date for the title;

last = length(xlab2);
dt = [ xlab1(1,1:8),' ',xlab2(1,1:5),' to ',xlab1(last,1:8),' ',xlab2(last,1:5)];


% find the max and mins for scaling and labeling the graphs

max3 = max(difference);
min3 = min(difference);

max_val = max([max1 max2 max3]);
min_val = min([min1 min2 min3]);
range = max_val - min_val; 

text(xar(5),max_val+0.13*range,dt,'FontSize',12,'FontWeight',...
  'bold','HorizontalAlignment','Center');
text(xar(5),max_val+0.21*range,'Difference Plot of Significant Wave Height',...
  'FontSize',14,'FontWeight','bold','HorizontalAlignment','Center');
text(xar(5),max_val+0.17*range,[char(info1.stnid),' ',char(info1.stname),' '],...
  'FontSize',14,'FontWeight','bold',...
  'Color','b','HorizontalAlignment','Right');
text(xar(5),max_val+0.17*range,[' ',char(info2.stnid),' ',char(info2.stname)],...
  'FontSize',14,'FontWeight','bold',...
 'Color',dg,'HorizontalAlignment','Left');

set(gca,'XLim',[xar(1) xar(last)]);
set(gca,'XTick',xar);
set(gca,'XTickLabel',xlab1);
set(gca,'Position',[ .12 .13 .8 .75 ],'units','normal');

sl=size(xlab1);
for i=1:sl(1)
  text(xar(i),min_val-.05*range,xlab2(i,1:5),'HorizontalAlignment','center');
end

ylim([min_val max_val*1.1]);
ylabel(pmlabel,'FontSize',12,'FontWeight','bold');
text(xar(5),min_val-.11*range,'Date (UTC)','HorizontalAlignment',...
  'center','FontSize',12,'FontWeight','bold');

% produce a legend with the appropriate colors

h = legend('Station 1','Station 2','Difference (1-2)',1);
ch = get(h,'Children');
%set(ch(1),'Color',[ 1 0 0 ],'LineWidth',4);
%set(ch(2),'Color',dg,'LineWidth',2);
%set(ch(3),'Color',[ 0 0 1 ],'LineWidth',2);

grid;

blue_box(0.01);

set(gcf, 'visible','off');
plt_png = [wkg,'plot_',pid,'.png'];
eval(['print -dpng -r75 ',plt_png]);

quit;