Commit bb4e183e authored by Pierre Cazenave's avatar Pierre Cazenave

Remove some broken scripts and some which were specific to given proejcts.

parent c71a1502
function [Plots]=do_ph_change_plot(plotOPTS,FVCOM,startIdx)
% Calculate the change in the given parameter (plotOPTS.var_plot) and plot
% accordingly.
m_mappath;
% For testing
% startIdx=48; % Leak start
% layerIdx=10; % Bottom layer
% For testing
% Check we have some of the required fields.
if ~isfield(FVCOM,plotOPTS.var_plot)
error('Need %s input to calculate change in %s.',plotOPTS.var_plot,plotOPTS.var_plot)
end
% Get the relevant time intervals
if isempty(plotOPTS.save_intervals)
plotOPTS.save_intervals=1:length(plotOPTS.Time_record);
dontSave=1;
else
dontSave=0;
end
% Build a colour palette which matches Jerry's ranges.
% Dark Red -> Red -> Amber -> Yellow -> Green -> Blue
nColours=200;
nColourIn=[1,nColours*0.15,nColours*0.6,nColours*0.75,nColours*0.9,nColours];
nColourOut=1:nColours; % Gives a nice continuous palette.
% DR R A Y G B
cRed=interp1(nColourIn, [0.62, 0.9, 1 , 1 , 0 , 0.46],nColourOut);
cGreen=interp1(nColourIn,[0 , 0 , 0.52 , 1 , 0.8, 0.63],nColourOut);
cBlue=interp1(nColourIn, [0.2 , 0.2, 0 , 0 , 0 , 0.83],nColourOut);
colourSpec=[cRed;cGreen;cBlue]';
if isfield(plotOPTS,'altColours') && plotOPTS.altColours==1
clear nColours nColourIn nColourOut colourSpec
% Build a colour palette which matches Jerry's ranges.
% Dark Red -> Red -> Amber -> Yellow -> Green -> Blue
nColours=200;
nColourIn=[1,nColours*0.1,nColours*0.2,nColours*0.3,nColours*0.4,nColours];
nColourOut=1:nColours; % Gives a nice continuous palette.
% DB B LB DG G LG
cRed=interp1(nColourIn, [0 , 0 , 0 , 0.04, 0 , 0.2 ],nColourOut);
cGreen=interp1(nColourIn,[0 , 0 , 0.2 , 0.51, 0.6, 0.76 ],nColourOut);
cBlue=interp1(nColourIn, [0.4 , 0.6 , 0.79, 0.78, 0.6 , 0],nColourOut);
colourSpec=flipud([cRed;cGreen;cBlue]');
end
dataToUse=single(FVCOM.(plotOPTS.var_plot));
% Get the background condition. Depth average if necessary.
if isfield(plotOPTS,'depth_average') && plotOPTS.depth_average
bgPH=squeeze(mean(dataToUse(:,:,startIdx),2));
else
bgPH=squeeze(dataToUse(:,plotOPTS.nz_plot,startIdx));
end
% Are we depth averaging?
if isfield(plotOPTS,'depth_average') && plotOPTS.depth_average
phData=squeeze(mean(dataToUse,2));
else
phData=squeeze(dataToUse(:,plotOPTS.nz_plot,:));
end
% Check if we're doing the cumulative difference. For this to work, you
% need to start your data import before the start of the leak.
% if isfield(plotOPTS,'summed_ph') && plotOPTS.summed_ph
% phDiff=cumsum(diff(phData,[],2));
% else
% % Otherwise, just do the difference from the start, ignoring all
% % previous steps i.e. for each successive time step, calculate the
% % difference between the current time step and the background
% % level.
% phDiff=phData-repmat(bgPH,1,size(phData,2));
% end
% Check if we're doing the cumulative difference
% if isfield(plotOPTS,'summed_ph') && plotOPTS.summed_ph
% phDiff=cumsum(phDiff,2);
% end
% Since we're not doing the cumulative difference, just do the difference
% relative to the background value
phDiff=phData-repmat(bgPH,1,size(phData,2));
figure(plotOPTS.figure); clf
m_proj('UTM','lon',[plotOPTS.range_lon],'lat',[plotOPTS.range_lat],'zon',plotOPTS.zone,'ell','grs80')
m_grid('box','fancy')
m_usercoast(plotOPTS.coastline_file,'Color','k','LineWidth',3);
[X,Y]=m_ll2xy(plotOPTS.mesh.lon,plotOPTS.mesh.lat,'clip','on');
saveInc=1;
for aa=1:length(plotOPTS.Time_record)
% plot map with change
hold on
Plots(1).handles=patch('Vertices',[X,Y],'Faces',plotOPTS.mesh.tri,...
'Cdata',phDiff(:,plotOPTS.save_intervals(aa)),...
'edgecolor','interp','facecolor','interp');
fprintf('Time step %i of %i\n',aa,length(plotOPTS.Time_record))
caxis(plotOPTS.clims)
colormap(colourSpec)
ch=colorbar;
set(ch,'FontSize',10);
ylabel(ch,'pH change');
% ylabel(ch,[plotOPTS.var_plot,' change'])
% check if mesh elements are required
if plotOPTS.do_mesh
% plot vertices
[X,Y]=m_ll2xy(plotOPTS.mesh.lon,plotOPTS.mesh.lat,'clip','on');
patch('Vertices',[X,Y],'Faces',plotOPTS.mesh.tri,...
'EdgeColor',[0.6 0.6 0.6],'FaceColor','none'); hold on
end
pause(plotOPTS.pause)
if plotOPTS.save_output % Are we even trying to save figures?
if saveInc<=length(plotOPTS.save_intervals) && plotOPTS.save_intervals(aa)==plotOPTS.save_intervals(saveInc) && ~dontSave
% Save output
fprintf('Saving figure... ')
set(findobj(gcf,'Type','text'),'FontSize',10)
%set(gcf,'PaperSize',fliplr(get(gcf,'PaperSize')))
set(gcf,'PaperPositionMode','auto');
set(gcf,'renderer','painters'); % for vector output in pdfs
print(gcf,'-dpdf','-r600',[plotOPTS.FVCOM_plot_dir,plotOPTS.var_plot,'/pdf/',plotOPTS.fig_name,'_layer=',num2str(plotOPTS.nz_plot),'_',plotOPTS.var_plot,'_change_',num2str(plotOPTS.save_intervals(aa)),'.pdf']); % pdf
%print(gcf,'-dpng','-r600',[plotOPTS.FVCOM_plot_dir,plotOPTS.var_plot,'/png/',plotOPTS.fig_name,'_layer=',num2str(plotOPTS.nz_plot),'_',plotOPTS.var_plot,'_change_',num2str(plotOPTS.save_intervals(aa)),'.png']); % png
saveInc=saveInc+1;
fprintf('done.\n')
end
end
if aa~=length(plotOPTS.Time_record)
delete(Plots(1).handles)
end
end
return
function [Plots]=do_ph_change_points_plot(plotOPTS,FVCOM,startIdx)
% Calculate the change in the given parameter (plotOPTS.var_plot) and plot
% accordingly.
m_mappath;
plotOPTS.figure=1;
plotOPTS.var_plot='ph';
plotOPTS.do_mesh=1; % Add the mesh?
startIdx=1; %:size(FVCOM.zeta,2);
plotOPTS.threshold_change=-0.2;
% near bottom, but not actually the bottom since it doesn't change at all
plotOPTS.nz_plot=1;
[nx,nz,ttot]=size(FVCOM.(plotOPTS.var_plot));
% Check we have some of the required fields.
if ~isfield(FVCOM,plotOPTS.var_plot)
error('Need %s input to calculate change in %s.',plotOPTS.var_plot,plotOPTS.var_plot)
end
% Get the relevant time intervals
if isfield(plotOPTS,'save_intervals') && isempty(plotOPTS.save_intervals)
plotOPTS.save_intervals=1:length(plotOPTS.Time_record);
dontSave=1;
else
dontSave=0;
end
% Get the time step (in seconds)
dt=round((plotOPTS.Time_record(2)-plotOPTS.Time_record(1))*24*60*60);
% Sigma layer fraction (nz*nz)
sigThickness=roundn(abs(diff(FVCOM.siglev,1,2)),-5); % roundn to even values out.
% Total depth (nz*ttot)
totalDepth=repmat(FVCOM.h,1,size(FVCOM.zeta,2))+FVCOM.zeta;
% Volume for every element for each time step
cellVolume=nan(nx,nz,ttot);
for ii=1:nz
for jj=1:ttot
% I'm pretty certain this needs to integrate for the time between
% time steps, otherwise halving the output time step would halve
% the exposed volume, which can't be right...
% However, if we know that the threshold has been triggered n
% times, we just need n times the volume, irrespective of time.
% Right?
cellVolume(:,ii,jj)=totalDepth(:,jj).*sigThickness(:,ii).*FVCOM.art1; %*dt;
end
end
% Calculate the difference between the current step and the previous step.
% phChange=diff(squeeze(FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,:)),[],2);
% Nope, get the background condition and compare against that. Depth
% average if necessary.
if isfield(plotOPTS,'depth_average') && plotOPTS.depth_average
bgPH=squeeze(mean(FVCOM.(plotOPTS.var_plot)(:,:,startIdx),2));
else
bgPH=squeeze(FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,startIdx));
end
if isfield(plotOPTS,'depth_average') && plotOPTS.depth_average
% Now, for each successive time step, calculate the difference between the
% current time step and the background level.
phMean=squeeze(mean(FVCOM.(plotOPTS.var_plot),2));
phDiff=phMean-repmat(bgPH,1,size(phMean,2));
else
phDiff=squeeze(FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,:))-repmat(bgPH,1,size(FVCOM.(plotOPTS.var_plot),3));
end
% Try a cumulative difference
phDiffCumulative=zeros(nx,ttot,'single');
for tt=1:ttot
if tt>1 && tt<ttot
currpH=FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,tt);
nextpH=FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,tt+1);
phDiffCumulative(:,tt)=phDiffCumulative(:,tt)+(currpH-nextpH);
elseif tt==1
currpH=FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,tt);
nextpH=FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,tt+1);
phDiffCumulative(:,tt)=currpH-nextpH;
end
end
[X,Y]=m_ll2xy(plotOPTS.mesh.lon,plotOPTS.mesh.lat,'clip','on');
foundXY=cell(1,length(plotOPTS.nz_plot));
countSites=zeros(nx,1);
for tt=1:ttot
% Find drops greater than the threshold
idx=find(phDiff(:,tt)<plotOPTS.threshold_change);
if ~isempty(idx)
countSites(idx)=countSites(idx)+1;
foundXY{tt}.xy=[X(idx),Y(idx)];
end
end
% Now we can use countSites to calculate the volume exposed to the change,
% although this ignores the effect of the tide (because we don't know which
% time step triggered the threshold condition test)
totalVolume=squeeze(cellVolume(:,plotOPTS.nz_plot,1)).*countSites;
% Visualise the count threshold
Plots(1).handles=patch('Vertices',[X,Y],'Faces',plotOPTS.mesh.tri,...
'Cdata',countSites,...
'edgecolor','interp','facecolor','interp');
if plotOPTS.do_mesh
% plot vertices
[X,Y]=m_ll2xy(plotOPTS.mesh.lon,plotOPTS.mesh.lat,'clip','on');
patch('Vertices',[X,Y],'Faces',plotOPTS.mesh.tri,...
'EdgeColor',[0.6 0.6 0.6],'FaceColor','none'); hold on
end
colorbar
figure;
patch('Vertices',[X,Y],'Faces',plotOPTS.mesh.tri,...
'Cdata',totalVolume/1e9,...
'edgecolor','interp','facecolor','interp');
if plotOPTS.do_mesh
% plot vertices
[X,Y]=m_ll2xy(plotOPTS.mesh.lon,plotOPTS.mesh.lat,'clip','on');
patch('Vertices',[X,Y],'Faces',plotOPTS.mesh.tri,...
'EdgeColor',[0.6 0.6 0.6],'FaceColor','none'); hold on
end
colorbar
figure;
patch('Vertices',[X,Y],'Faces',plotOPTS.mesh.tri,...
'Cdata',squeeze(FVCOM.ph(:,1,2)),...
'edgecolor','interp','facecolor','interp');
colorbar
% figure(1)
% subplot(2,1,1)
% plot(Time_record-min(Time_record),squeeze(FVCOM.ph(1316,plotOPTS.nz_plot,:)))
% subplot(2,1,2)
% plot(Time_record(1:end-1)-min(Time_record),phChange(1316,:),'r')
% hold on
% plot(Time_record(idxAtLeak)-min(Time_record),phChange(1316,idxAtLeak),'g.')
% % Add the threshold
% plot([min(Time_record)-min(Time_record),max(Time_record)-min(Time_record)],[plotOPTS.threshold_change,plotOPTS.threshold_change],'k--')
function [Plots]=do_ph_max_change_plot(plotOPTS,FVCOM,startIdx)
% Calculate the change in pH and plot accordingly.
m_mappath;
% Check we have some of the required fields.
if ~isfield(FVCOM,plotOPTS.var_plot)
error('Need %s input to calculate change in %s.',plotOPTS.var_plot,plotOPTS.var_plot)
end
% Build a colour palette which matches Jerry's ranges.
% Dark Red -> Red -> Amber -> Yellow -> Green -> Blue
nColours=200;
nColourIn=[1,nColours*0.15,nColours*0.6,nColours*0.75,nColours*0.9,nColours];
nColourOut=1:nColours; % Gives a nice continuous palette.
% DR R A Y G B
cRed=interp1(nColourIn, [0.62, 0.9, 1 , 1 , 0 , 0.46],nColourOut);
cGreen=interp1(nColourIn,[0 , 0 , 0.52 , 1 , 0.8, 0.63],nColourOut);
cBlue=interp1(nColourIn, [0.2 , 0.2, 0 , 0 , 0 , 0.83],nColourOut);
colourSpec=flipud([cRed;cGreen;cBlue]'); % flip since we have negative scale
if isfield(plotOPTS,'altColours') && plotOPTS.altColours==1
clear nColours nColourIn nColourOut colourSpec
% Build a colour palette which matches Jerry's ranges.
% Dark Red -> Red -> Amber -> Yellow -> Green -> Blue
nColours=200;
nColourIn=[1,nColours*0.1,nColours*0.2,nColours*0.3,nColours*0.4,nColours];
nColourOut=1:nColours; % Gives a nice continuous palette.
% DB B LB DG G LG
cRed=interp1(nColourIn, [0 , 0 , 0 , 0.04, 0 , 0.2 ],nColourOut);
cGreen=interp1(nColourIn,[0 , 0 , 0.2 , 0.51, 0.6, 0.76 ],nColourOut);
cBlue=interp1(nColourIn, [0.4 , 0.6 , 0.79, 0.78, 0.6 , 0],nColourOut);
colourSpec=[cRed;cGreen;cBlue]';
end
dataToUse=single(FVCOM.(plotOPTS.var_plot));
% Get the background condition. Depth average if necessary.
if isfield(plotOPTS,'depth_average') && plotOPTS.depth_average
bgPH=squeeze(mean(dataToUse(:,:,startIdx),2));
else
bgPH=squeeze(dataToUse(:,plotOPTS.nz_plot,startIdx));
end
% Are we depth averaging?
if isfield(plotOPTS,'depth_average') && plotOPTS.depth_average
phData=squeeze(mean(dataToUse,2));
else
phData=squeeze(dataToUse(:,plotOPTS.nz_plot,:));
end
phDiff=phData-repmat(bgPH,1,size(phData,2));
% Are we depth averaging?
% if isfield(plotOPTS,'depth_average') && plotOPTS.depth_average
% % Calculate the difference for each time step from the initial
% % conditions.
% phDiff=squeeze(mean(FVCOM.(plotOPTS.var_plot),2))-repmat(bgPH,1,length(plotOPTS.Time_record));
% else
% phDiff=squeeze(FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,:))-repmat(bgPH,1,size(FVCOM.(plotOPTS.var_plot),3));
% end
% Not doing cumulative differences any more
% Check if we're doing the cumulative difference
% if isfield(plotOPTS,'summed_ph') && plotOPTS.summed_ph
% phDiff=cumsum(phDiff,2);
% end
% Do the maximum change from initial conditions across all time steps (i.e.
% not the cumulative change). abs() the difference to find the greatest
% change, either increase or decrease.
% We're also averaging here over the number of timesteps so that we can fix
% the colour palette to something sensible.
phDiffMax=min((phDiff),[],2); %./(size(plotOPTS.Time_record,2)*(plotOPTS.Time_record(2)-plotOPTS.Time_record(1)));
figure(plotOPTS.figure); clf
m_proj('UTM','lon',[plotOPTS.range_lon],'lat',[plotOPTS.range_lat],'zon',plotOPTS.zone,'ell','grs80')
m_grid('box','fancy')
m_usercoast(plotOPTS.coastline_file,'Color','k','LineWidth',3);
[X,Y]=m_ll2xy(plotOPTS.mesh.lon,plotOPTS.mesh.lat,'clip','on');
% plot map with pH change
hold on
Plots(1).handles=patch('Vertices',[X,Y],'Faces',plotOPTS.mesh.tri,...
'Cdata',phDiffMax,'edgecolor','interp','facecolor','interp');
caxis(plotOPTS.clims)
colormap(flipud(colourSpec))
ch=colorbar;
set(ch,'FontSize',10);
ylabel(ch,'pH change');
% ylabel(ch,[plotOPTS.var_plot,' change'])
% check if mesh elements are required
if plotOPTS.do_mesh
% plot vertices
[X,Y]=m_ll2xy(plotOPTS.mesh.lon,plotOPTS.mesh.lat,'clip','on');
patch('Vertices',[X,Y],'Faces',plotOPTS.mesh.tri,...
'EdgeColor',[0.6 0.6 0.6],'FaceColor','none'); hold on
end
if plotOPTS.save_output % Are we even trying to save figures?
% Save output
fprintf('Saving figure... ')
set(findobj(gcf,'Type','text'),'FontSize',10)
%set(gcf,'PaperSize',fliplr(get(gcf,'PaperSize')))
set(gcf,'PaperPositionMode','auto');
set(gcf,'renderer','painters'); % for vector output in pdfs
print(gcf,'-dpdf','-r600',[plotOPTS.FVCOM_plot_dir,plotOPTS.var_plot,'/pdf/',plotOPTS.fig_name,'_layer=',num2str(plotOPTS.nz_plot),'_',plotOPTS.var_plot,'_max_change.pdf']); % pdf
%print(gcf,'-dpng','-r600',[plotOPTS.FVCOM_plot_dir,plotOPTS.var_plot,'/png/',plotOPTS.fig_name,'_layer=',num2str(plotOPTS.nz_plot),'_',plotOPTS.var_plot,'_max_change.png']); % png
fprintf('done.\n')
end
return
function [Plots]=do_ph_vertical_profile(plotOPTS,FVCOM,transectPoints)
% For a supplied transect, plot the vertical profile through the water
% column supplied in FVCOM.var_plot.
% Get the relevant time intervals
if isempty(plotOPTS.save_intervals)
plotOPTS.save_intervals=1:length(plotOPTS.Time_record);
dontSave=1;
else
dontSave=0;
end
if isfield(plotOPTS,'altColours') && plotOPTS.altColours==1
colourSpec=flipud(jet);
else
clear nColours nColourIn nColourOut colourSpec
% Build a colour palette which matches Jerry's ranges.
% Dark Red -> Red -> Amber -> Yellow -> Green -> Blue
nColours=200;
nColourIn=[1,nColours*0.15,nColours*0.6,nColours*0.75,nColours*0.9,nColours];
nColourOut=1:nColours; % Gives a nice continuous palette.
% DR R A Y G B
cRed=interp1(nColourIn, [0.62, 0.9, 1 , 1 , 0 , 0.46],nColourOut);
cGreen=interp1(nColourIn,[0 , 0 , 0.52 , 1 , 0.8, 0.63],nColourOut);
cBlue=interp1(nColourIn, [0.2 , 0.2, 0 , 0 , 0 , 0.83],nColourOut);
colourSpec=[cRed;cGreen;cBlue]';
end
dataToUse=single(FVCOM.(plotOPTS.var_plot));
saveInc=1;
for aa=1:length(plotOPTS.Time_record)
verticalProfile=squeeze(dataToUse(transectPoints.trn_nodes.idx,:,plotOPTS.save_intervals(aa)));
% xValues=repmat(transectPoints.trn_dis,1,size(verticalProfile,2));
% Not sure about those distance values...
xNorm=transectPoints.trn_nodes.x-min(transectPoints.trn_nodes.x);
yNorm=transectPoints.trn_nodes.y-min(transectPoints.trn_nodes.y);
xDist=sqrt(xNorm.^2+yNorm.^2);
xValues=repmat(xDist,1,size(verticalProfile,2));
zeta=squeeze(FVCOM.zeta(transectPoints.trn_nodes.idx,plotOPTS.save_intervals(aa)));
waterDepth=FVCOM.h(transectPoints.trn_nodes.idx);
zValues=(zeta+waterDepth)*FVCOM.siglay(1,:);
% plot profile
fprintf('Time step %i of %i\n',aa,length(plotOPTS.Time_record))
figure(plotOPTS.figure); clf
hold on
Plots(1).handles=contourf(...
xValues/1000,... % km
zValues,...
fliplr(verticalProfile),... % get vertical right way up
200,'edgecolor','none');
% ylim([-max(FVCOM.h(:)) -0.5])
xlabel('Distance along transect (km)')
ylabel('Depth (m)')
colormap(colourSpec)
ch=colorbar;
set(ch,'FontSize',10);
ylabel(ch,'pH')
if ~isempty(plotOPTS.vlims)
caxis(plotOPTS.vlims);
% Get sensible tick formatting (if necessary)
if plotOPTS.vlims(2)-plotOPTS.vlims(1)<1e-4
xx=0:length(plotOPTS.vlims)-1;
yy=plotOPTS.vlims;
xxi=0:1/6:1; % six ticks
cticks=interp1(xx,yy,xxi);
set(ch,'YTick',double(cticks))
set(ch,'yticklabel',cticks)
end
end
pause(plotOPTS.pause)
if plotOPTS.save_output % Are we even trying to save figures?
if saveInc<=length(plotOPTS.save_intervals) && plotOPTS.save_intervals(aa)==plotOPTS.save_intervals(saveInc) && ~dontSave
% Save output
fprintf('Saving figure... ')
set(findobj(gcf,'Type','text'),'FontSize',10)
%set(gcf,'PaperSize',fliplr(get(gcf,'PaperSize')))
set(gcf,'PaperPositionMode','auto');
set(gcf,'renderer','painters'); % for vector output in pdfs
print(gcf,'-dpdf','-r600',[plotOPTS.FVCOM_plot_dir,plotOPTS.var_plot,'/pdf/',plotOPTS.fig_name,'_vertical_profile_',plotOPTS.var_plot,'_',num2str(plotOPTS.save_intervals(aa)),'.pdf']); % pdf
%print(gcf,'-dpng','-r600',[plotOPTS.FVCOM_plot_dir,plotOPTS.var_plot,'/png/',plotOPTS.fig_name,'_layer=',num2str(plotOPTS.nz_plot),'_',plotOPTS.var_plot,'_change_',num2str(plotOPTS.save_intervals(aa)),'.png']); % png
saveInc=saveInc+1;
fprintf('done.\n')
end
end
if aa~=length(plotOPTS.Time_record)
clf
end
end
return
\ No newline at end of file
function [Plots,totalVol]=do_volume(plotOPTS,FVCOM,startIdx,thresholdValue)
% Calculate a volume for a given time step above a threshold value.
m_mappath;
[nx,nz,ttot]=size(FVCOM.(plotOPTS.var_plot));
% Get the time step (in seconds)
dt=round((plotOPTS.Time_record(2)-plotOPTS.Time_record(1))*24*60*60);
% Sigma layer fraction (nz*nz)
sigThickness=roundn(abs(diff(FVCOM.siglev,1,2)),-5); % roundn to even values out.
% Total depth (nz*ttot)
totalDepth=repmat(FVCOM.h,1,size(FVCOM.zeta,2))+FVCOM.zeta;
% Volume for every element for each time step
cellVolume=nan(nx,nz,ttot);
for ii=1:nz
for jj=1:ttot
cellVolume(:,ii,jj)=totalDepth(:,jj).*sigThickness(:,ii).*FVCOM.art1*dt;
end
end
if plotOPTS.nz_plot==0
% Depth averaging the results
cellVolume=squeeze(sum(cellVolume,2));
dataArray=squeeze(mean(FVCOM.(plotOPTS.var_plot)(:,:,:),2));
else
dataArray=FVCOM.(plotOPTS.var_plot)(:,:,:);
end
%%
% Calculate volume of elements matching threshold condition
totalVolTemp=0;
colourSpec=hsv(length(plotOPTS.nz_plot));
plotSymbols={'.','o','x','+','^','*','p','h'};
for tt=1:size(startIdx,2)
foundXY=cell(1,length(plotOPTS.nz_plot));
figure(plotOPTS.figure); clf
m_proj('UTM','lon',[plotOPTS.range_lon],'lat',[plotOPTS.range_lat],'zon',plotOPTS.zone,'ell','grs80')
m_grid('box','fancy')
m_usercoast(plotOPTS.coastline_file,'Color','k','LineWidth',3);
[X,Y]=m_ll2xy(plotOPTS.mesh.lon,plotOPTS.mesh.lat,'clip','on');
fprintf('Time step %i of %i\n',startIdx(tt),length(plotOPTS.Time_record))
% Find all values of pH greater than the thresholdValue for the layers
% specified in .nz_plot only.
for ss=1:size(plotOPTS.nz_plot,2)
if plotOPTS.nz_plot==0
% Depth averaged, so easy.
idx=thresholdTest(dataArray(:,startIdx(tt)),plotOPTS.volume_thresh,thresholdValue);
% Do we have any locations where the condition has been met?
if ~isempty(idx)
totalVolTemp=totalVolTemp+sum(sum(cellVolume(idx,startIdx(tt))));
end
else
% Loop through the layers
for qq=1:size(plotOPTS.nz_plot,2)
idx=thresholdTest(dataArray(:,plotOPTS.nz_plot(qq),startIdx(tt)),plotOPTS.volume_thresh,thresholdValue);
% Again, look for areas where we have met the threshold
% condition. Here, we need to save the xy coordinates for
% each find with the layer number.
if ~isempty(idx)
totalVolTemp=totalVolTemp+sum(sum(cellVolume(idx,plotOPTS.nz_plot(qq),startIdx(tt))));
% This may skip the odd layer...
foundXY{qq}.xy=[X(idx),Y(idx),repmat(plotOPTS.nz_plot(qq),length(idx),1)];
end
end
end
end
if tt==1;
totalVol(tt)=totalVolTemp;
else
totalVol(tt)=totalVol(tt-1)+totalVolTemp;
end
% Do a figure identifying the locations above the threshold for the
% surface sigma layer.
if plotOPTS.nz_plot==0
Plots(plotOPTS.figure).handles=patch('Vertices',[X,Y],'Faces',...
plotOPTS.mesh.tri,'Cdata',squeeze(cellVolume(:,tt)),...
'edgecolor','interp','facecolor','interp');
else
Plots(plotOPTS.figure).handles=patch('Vertices',[X,Y],'Faces',...
plotOPTS.mesh.tri,'Cdata',squeeze(cellVolume(:,1,tt)),...
'edgecolor','interp','facecolor','interp');
end
hold on
% TODO: For multiple layers, do unique indices for each layer so they
% can be plotted separately here (a la do_vector_plot.m).
if ~isempty(idx)
for kk=1:size(plotOPTS.nz_plot,2);
plot(foundXY{kk}.xy(:,1),foundXY{kk}.xy(:,2),plotSymbols{mod(kk,length(plotSymbols)+1)},'Color',colourSpec(kk,:))
end
end
if plotOPTS.do_mesh
% plot vertices
[X,Y]=m_ll2xy(plotOPTS.mesh.lon,plotOPTS.mesh.lat,'clip','on');
patch('Vertices',[X,Y],'Faces',plotOPTS.mesh.tri,...
'EdgeColor',[0.6 0.6 0.6],'FaceColor','none'); hold on
end
% Some useful text
[textX,textY]=m_ll2xy(-4.45,50.13);
text(textX,textY,sprintf('Total volume:\t\t%.2fkm^{3}\nVolume this time step:\t%.2fkm^{3}\n',totalVol(tt)/1e9,totalVolTemp/1e9))
pause(plotOPTS.pause)
caxis(plotOPTS.clims)
colorbar
set(get(colorbar,'YLabel'),'String','Volume (m^{3})')
if tt~=size(startIdx,2)
delete(Plots(plotOPTS.figure).handles)
end
end
function idx=thresholdTest(dataArray,thresholdType,thresholdValue)
% Abstract some of the complexity into a function.
% Usage: thresholdTest(inputData,thresholdValue)
switch thresholdType
case -1
idx=find(dataArray<thresholdValue);
case 1
idx=find(dataArray>thresholdValue);
case 0
idx=find(dataArray==thresholdValue);
otherwise
error('Unrecognised value for ''plotOPTS.volume_thresh''.')
end
return
\ No newline at end of file
function [totalVolume] = do_volume_change(plotOPTS,FVCOM)
% DO_VOLUME_CHANGE Calculate volume of water which experiences a change in
% pH beyond some defined threshold.
% Make a 3D array of the volumes to get to total volume which experiences
% the change in pH.
[nx,nz,ttot]=size(FVCOM.(plotOPTS.var_plot));
% Get the time step (in seconds)
dt=round((plotOPTS.Time_record(2)-plotOPTS.Time_record(1))*24*60*60);
% Make sure we're actually able to sample at the rate requested.
if plotOPTS.change_type~=0 && dt/(60*60)>plotOPTS.change_type
error('Output file sampling frequency is coarser than specified time sampling.')
end
% Sigma layer fraction (nz*nz)
sigThickness=roundn(abs(diff(FVCOM.siglev,1,2)),-5); % roundn to even values out.
% Total depth (nz*ttot)
totalDepth=repmat(FVCOM.h,1,size(FVCOM.zeta,2))+FVCOM.zeta;
% Volume for every element for each time step
cellVolume=nan(nx,nz,ttot);
for ii=1:nz
for jj=1:ttot
cellVolume(:,ii,jj)=totalDepth(:,jj).*sigThickness(:,ii).*FVCOM.art1.*dt;
end
end
if plotOPTS.change_type==0 || plotOPTS.change_type==dt/(60*60)
% Instantaneous change (i.e. a change between outputs beyond a given value.
phChange=diff(FVCOM.DYE,1,3);
% Find the elements with the change beyond the threshold.
idx=phChange<plotOPTS.threshold_change;
% Using those indices, find the volume which is affected.
totalVolume=sum(cellVolume(idx));
else
% The more complicated n-hour change.
% Use the time period over which we're interested in to get an index
% jump value.
dtJump=(plotOPTS.change_type*60*60)/dt;
if dtJump-round(dtJump)~=0
error('Output sampling is not compatible with time period of change.')
end
phChange=nan(nx,nz,ttot-dtJump);
for tt=1:length(plotOPTS.Time_record)
if tt<=length(plotOPTS.Time_record)-dtJump
phChange(:,:,tt)=(FVCOM.DYE(:,:,tt+dtJump)-FVCOM.DYE(:,:,tt));
end
end
% Find the change values that happen for at least n-hours.
% Probably best to look away now if you're interested in clean quick
% code.
totalVolume=0;
for ii=1:nx
for jj=1:nz
totalVolumeCurrent=get_runs(plotOPTS,squeeze(phChange(ii,jj,:)),squeeze(cellVolume(ii,jj,:)),plotOPTS.change_type,plotOPTS.threshold_change);
totalVolume=totalVolume+totalVolumeCurrent;
end
end
end
end
\ No newline at end of file
function [summedVolume]=do_volume_change2(plotOPTS,FVCOM,startIdx)
% Calculate the volume exposed to a given change in pH.
m_mappath;