do_ph_max_change_plot.m 4.69 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
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