ENERGY SPECTRUM MONTHLY PLOT
Energy spectrum data - derived from a fast fourier transform of approximately 30 minutes of wave data - is plotted with frequency on the x axis and time on the y axis. The elevation of a "mountain" corresponds to the energy density, with the scale to the right of the plot. Every fourth data run is plotted with one of four colors.
This plot is useful for observing the evolution of wave energy over all the measured frequencies during the month.
Matlab code for generating the plot (no longer used)
%-------------------------------------------------------------------------------
% Name:
% .dbase/web_programs/plotpro/matlab_scripts/mountain.m
% Description:
% Creates the monthly spectral energy plot (mountains).
% A few lines are pre-pended to this script.
% Called by:
% .f90/plot_utils.make_plot
%-------------------------------------------------------------------------------
path(path,'/project/dbase/web_programs/plotpro/matlab_scripts');
wkg = wkg_path;
% LOAD DATA, GIVE SHORT NAMES FOR LOADED VARIABLES
% NOTE: x = frequency data, y = energy density data
eval(['freq=load(''',wkg,'freq.',pid,'01'');']);
eval(['ener=load(''',wkg,'ener.',pid,'01'');']);
x = freq(1,:);
y = ener;
dt = read_date_file([pid,'01'],wkg);
tim = (dt.day-1)*24 + (dt.hour + dt.min/60); % in hours
yearchar = num2str(dt.year(1));
monthnum = dt.month(1);
clear freq ener date
% CHECK TO SEE IF THERE IS ANY DATA
if ~ ( any(x>0) & any(y>0) )
clear all;
quit;
end
% SET FIGURE WINDOW
figure('units','normal','paperunits','normal','position',...
[.1,.1,.5,.8],'paperposition',[.05,.05,.9,.93],'visible','off');
sz = size(y);
map = [ 0 0 0; 1 0 0; 0 1 0; 0 0 1; .5 .5 .5];
colormap(map);
cd = colormap;
% TAG FIRST RUN OF EACH DAY
pday = 0;
for i=1:sz(1)
day = floor(tim(i)/24) + 1;
if day ~= pday
nday(i) = 1;
pday = day;
else
nday(i) = 0;
end
end
% SCALE Y VALUES
yscale = 7;
ym = max(max(y));
% Fixed scale: set ym=50.;
%ym=50.;
y = yscale*24*y./ym;
logx = log(x);
% PLOT DATA
eps = 5;
prvtime = 790;
j = 0;
for i=sz(1):-1:1
if ( (prvtime - tim(i) > eps) | nday(i) == 1)
patch('xdata',logx,'ydata',tim(i)+y(i,:),'edgecolor','none','facecolor',cd(1+rem(j,5),:));
j = j + 1;
prvtime = tim(i);
end
end
yaxisheight = .76;
set(gca,'position',[0.1 0.1 0.7 yaxisheight]);
% GET STATION INFO
info = read_info_file([pid,'01'],wkg);
% IDENTIFY IF SURGE BASIN OR STANDARD PROCESSING
surge = 0;
basin = 0;
standard = 0;
if ~isempty( findstr(info.stream_name,'SURGE') )
surge = 1;
elseif ~isempty( findstr(info.stream_name,'BASIN') )
basin = 1;
else
standard = 1;
end
% WRITE STATION INFO ON PLOT
if ( surge == 1 )
xtext=log(1/128);
elseif ( basin == 1 )
xtext=log(1/22);
else
xtext=log(.10);
end
ytext= 860;
txtspace = 24;
text(xtext,ytext+3*txtspace,[char(info.stnid),' ',char(info.stname)],...
'horizontalalignment','center','fontsize',14);
stream_name = ['(',char(info.stream_name),')'];
if ( length(stream_name) > 5 )
text(xtext, ytext+2*txtspace, stream_name, 'horizontalalignment',...
'center','fontsize',12);
end
mname = get_monthname(monthnum,'long');
text(xtext,ytext + -0.3*txtspace,[char(mname),' ',yearchar],'horizontalalignment',...
'center','fontsize',14);
text(xtext,ytext+0.7*txtspace,'ENERGY SPECTRUM','horizontalalignment','center',...
'fontsize',14);
% MAKE Y AXIS LABELS
ylim = 824;
set(gca,'ylim',[0,ylim],'ytick',[0:24:720],'yticklabel',num2str([1:31]'));
ylabel('Day of Month (UTC)','fontsize',12);
% MAKE X AXIS LABELS
numticks = 5;
if ( surge == 1 )
lxlim = 1/1024; % (Hz)
uxlim = 1/16; % (Hz)
f_sig_digs = 2;
p_sig_digs = 4;
f_labels = 0.20;
elseif ( basin == 1 )
lxlim = 1/128; % (Hz)
uxlim = 1/4; % (Hz)
f_sig_digs = 2;
p_sig_digs = 3;
f_labels = 0.20;
else
lxlim = 1/25; % (Hz)
uxlim = 1/4; % (Hz)
f_sig_digs = 2;
p_sig_digs = 2;
f_labels = 0.07;
end
% note: xtl = "x tick label"
diff = (log(uxlim) - log(lxlim))/(numticks-1);
xt = log(lxlim) + (0:(numticks-1))*diff;
xt(numticks) = log(uxlim);
fxtl = num2str(exp(xt'),f_sig_digs);
pxtl = num2str(1./exp(xt'),p_sig_digs);
set(gca,'xlim',[log(lxlim),log(uxlim)],'xtick',xt,'xticklabel',fxtl);
hold on;
plot([log(lxlim) log(lxlim)],[0 840],'k-');
box on;
text(xt,-35*ones(1,numticks),pxtl,'horizontalalignment','center');
text(log(uxlim)+f_labels,-15,'Frequency (Hz)','horizontalalignment','left',...
'fontsize',10);
text(log(uxlim)+f_labels,-35,'Period (s)','horizontalalignment','left',...
'fontsize',10);
% PUT SCALE ON PLOT
if ym < .1
rndval = .01;
elseif ym < 1
rndval = .1;
elseif ym < 10
rndval = 1;
else
rndval = 10;
end
rnd = rndval*round(ym/rndval)
yscaleheight = (yaxisheight/ylim) * (rnd/ym)*yscale*24;
sh=subplot('position',[.87,.60,.01,yscaleheight]);
hold on;
%draw the scale line (3 units long)
ytick = [0,rnd/4,rnd/2,3*rnd/4,rnd];
yticklabs = num2str(ytick',4);
plot([0,0],[0,rnd]);
hold on;
set(gca,'ylim',[0,rnd],'ytick',ytick,'yticklabel',yticklabs,'xticklabel',[]);
text(+9,rnd/2,'Energy Density (m^2/Hz)','rotation',90,'horizontalalignment','center','fontsize',12);
% PUT BLUE BOX AROUND PLOT AND ADD WEB ADDRESS
blue_box(0.01);
[pr_axes,ax_orig]=page_relative_axes;
axes(pr_axes);
text(.1,.035,'(Colors delineate data records)','FontSize',11);
axes(ax_orig);
% PRINT THE PLOT
set(gcf, 'visible','off');
png_name = [wkg,'plot_',pid,'.png'];
print('-dpng','-r75',png_name);
quit
