Add variable baseline functionality to anomaly plot function

[JakeWeis]

The anomaly plot function could be much more versatile by allowing the user to specify a variable baseline. This is useful for visualising differences between the annual trend of a given parameter and its climatological mean. This is an example plot:

Implementing this into the code is pretty straightforward:

1. Line 81

Assert that the base is either a scalar (constant baseline) or a vector the same length as x and y (variable baseline):

assert(isscalar(base) | numel(base)==numel(y),'Input error: Base value must be a scalar or a vector of the same length as x and y.')

2. Data Manipulation section (line 103)

If base is a scalar, convert it into a column vector the length of y (before "Columnate inputs to ensure...":

% If base is a scalar convert it into a column vector the length of y:
if isscalar(base)
    base = repmat(base,numel(y),1);
end

Columnate base just like x and y (only relevant if base is not a scalar):

base = base(:);

Find NaNs both in y and base so filling will work (and remove them from base as well):

% If y or base contains nans, ignore them so filling will work:
ind = (isfinite(y) & isfinite(base));
...
base = base(ind);

The intersections subfunction now simplifies as follows:

[xc,yc] = intersections(x,y,x,base);

Add zero crossings to base:

% Add zero crossings to the input dataset and sort them into the proper order: 
...
base = [base;yc];
...
base = base(ind);

Splitting the data into a top and bottom dataset changes as follows:

% Start thinking about this as two separate datasets which share baseline values where they meet: 
ytop = y; 
ytop(ytop<base) = base(ytop<base); 
ybot = y; 
ybot(ybot>base) = base(ybot>base); 

3. Plotting

Instead of using area we now need to use fill, the output remains the same:

% Plot the top half: 
htop = fill([x;flipud(x)],[ytop;flipud(base)],topcolor);

% Plot the bottom half: 
hbot = fill([x;flipud(x)],[ybot;flipud(base)],bottomcolor);

VOILA! This will allow the user to provide either a constant or a variable baseline and the code will do the rest:

Here is a copy of the adjusted code:
anomaly_JW.m.zip

I hope this is helpful.

Cheers,
Jake

[chadagreene]

Hi Jake,

This is a fantastic idea, and your implementation is brilliant. Thanks for describing the changes so meticulously.

Is there a reason you got rid of the threshold option? That's a bit of functionality I'd like to keep, because it's a somewhat common way of displaying ENSO anomalies

[JakeWeis]

Oops, that was not intentional. I was using the anomaly code published on File Exchange, which does not seem to have the threshold option. I'll rewrite it using the CDT anomaly function, which seems to be the most recent one.

[JakeWeis]

Alright, here we go:

The following edits were applied to the anomaly.m file in this repository:

Line 70

assert(numel(thresh)<=2 | numel(thresh)==numel(y),'Input error: The threshold must either be one or two scalars or the length of y.')

Data manipulation
The changes in this section come down to the following: I start off by converting the threshold input into a top and a bottom threshold column vector the size of y. In the following, this will allow us to handle the threshold the same for any input (single/double/variable threshold).

% Convert thresh into a top and a bottom column vector the length of y
if numel(thresh) == 1
    thresht = repmat(thresh,numel(y),1);
    threshb = repmat(thresh,numel(y),1);
elseif numel(thresh) == 2
    thresht = repmat(max(thresh),numel(y),1);
    threshb = repmat(min(thresh),numel(y),1);
else
    thresht = thresh(:);
    threshb = thresh(:);
end

This section is now redundant:

% if isscalar(thresh)
%    % Add a horizontal line at the end: 
%    x_archive = [x;NaN;min(x);max(x)]; 
%    y_archive = [y;NaN;thresh;thresh];
%    
%    % Ensure thresh has a lower and upper value: 
%    thresh = [thresh thresh]; 
% end
% thresh = sort(thresh); 

Check for and remove NaNs both in y and thresh:

% If y contains nans, ignore them so filling will work: 
ind = (isfinite(y) & isfinite(thresht) & isfinite(threshb)); 
...
thresht = thresht(ind); 
threshb = threshb(ind); 

The intersections functions change as follows (and will work for any threshold input):

[xct,yct] = intersections(x,y,x,thresht); % intersections is a subfunction by Douglas Schwarz, included below.
[xcb,ycb] = intersections(x,y,x,threshb); % intersections is a subfunction by Douglas Schwarz, included below.

Adding zero crossings to the both thresh vectors, and sorting them with xb/xt

% Add zero crossings to the input dataset and sort them into the proper order: 
...
thresht = [thresht;yct];
threshb = [threshb;ycb];
...
threshb = threshb(ind); % sorts threshb with xb
...
thresht = thresht(ind); % sorts thresht with xt

Separating y into a top and a bottom dataset changes as follows:

% Start thinking about this as two separate datasets which share threshline values where they meet: 
yb(yb>threshb) = threshb(yb>threshb); 
yt(yt<thresht) = thresht(yt<thresht); 

Plotting
As before, instead of using area we now need to use fill:

htop = fill([xt;flipud(xt)],[yt;flipud(thresht)],topcolor);
hbot = fill([xb;flipud(xb)],[yb;flipud(threshb)],bottomcolor);

This should be it :)

Code: anomaly_varthresh.m.zip
Sample data: sampledata.mat.zip