How to create a "skill-bias diagram" (meteorology)?

Question

In my research area (meteorology) graphs within graphs are commonly produced.

more information about it can be found here.

Each of those lines joints up data points that have:

1. An x-value, between 0 and 1 (values greater than 1 should not be represented in the graph).
2. A y-value, between 0 and 1.
3. A PSS value, between 1 and -1.
4. A Frequency Bias value, ranging from 0 to +∞, but values higher than 4 are not displayed.
5. A False Alarm Ratio (FAR) value, ranging from 0.0 to 0.9. The False Alarm Ratio value is held constant at a particular value for each data point on any given line.

EDIT: To make things really concrete, I've drawn a pink dot on the graph. That dot represents a data point for which x=0.81, y=0.61, PSS=-0.2, B=3.05, FAR=0.8.

I am trying to reproduce something similar in MATLAB. Googling turned up a lot of answers like this, which feature inset figures rather than what I'm looking for.

I have the data organized in a 3D array, where each page refers to a different level of False Alarm Ratio. The page with a FAR of 0.8 (data here) starts out like this

Then there are other pages on the 3D array devoted to FARs of 0.7, 0.6, and so on.

Questions
1. Is it even possible to create such an graph in MATLAB?
2. If so, what function should I use, and what approach should I take? EDIT: I have working code (below) that creates a somewhat similar figure using the linear plot function, but the documentation for this function does not indicate any way to insert a graph inside another graph. I am not sure how helpful this code is, but have inserted it in response to the downvoter.

H = [0:0.01:1];    
figure; hold on
fill([0 1 1],[0 0 1],[0 0.2 0.4]) % Deep blue
fill([0 1 0],[0 1 1],[0.4 0 0]) % Purple    
low_colours = {[0 0.501 1],[0 0.8 0.4], [0.4 0.8 0], [0.8 0.8 0]};
high_colours = {[0.6 0 0],[0.8 0 0], [1 0.5019 0], [0.988 0.827 0.196]};    
colour_counter = 0;

for ii = -0.8:0.2:0
    colour_counter = colour_counter + 1;
    if colour_counter < 5
        colour_now = low_colours{colour_counter};
    end
    ORSS = ones(1,size(H,2))*ii;
    F = (H .* (1-ORSS)) ./ ((1-2.*H) .* ORSS + 1);    
    plot(F,H)
    fill(F,H,colour_now);        
end

colour_counter = 0;

for ii = 0.8:-0.2:0
    colour_counter = colour_counter + 1;
    if colour_counter < 5
        colour_now = high_colours{colour_counter};
    end
    ORSS = ones(1,size(H,2))*ii;
    F = (H .* (1-ORSS)) ./ ((1-2.*H) .* ORSS + 1);
    plot(F,H)
    fill(F,H,colour_now);    
end

Answer 1

I think I got what you want, but before you go to the code below, notice the following:

1. I didn't need any of your functions in the link (and I have no idea what they do).
2. I also don't really use the x and y columns in the data, they are redundant coordinates to the PSS and B.
3. I concat all the 'pages' in your data to one long table (FAR below) with 5 columns (FAR,x,y,PSS,FB).

If you take closer look at the data you see that some areas that supposed to be colored in the graph has no representation in it (i.e. no values). So in order to interpolate the color to there we need to add the corners:

FAR{end+1,:} = [0.8 0 0 0 4];
FAR{end+1,:} = [0.9 0 0 -0.66 3.33];
FAR{end+1,:} = [1 0 0 0 0];
FAR{end+1,:} = [1 0 0 -1 3];

Next, the process has 2 parts. First we make a matrix for each variable, that ordered in columns by the corresponding FAR value, so for instance, in the PSS matrix the first column is all PSS values where FAR is 0, the second column is all PSS values where FAR is 0.1, and so on. We make such matrices for FAR(F), PSS and FreqBias(B), and we initialize them with NaNs so we can have columns with different number of values:

F = nan(max(histcounts(FAR.FAR,10)),10);
PSS = F;
B = F;
c = 1;
f = unique(FAR.FAR).';
for k = f
    valid = FAR.FAR==k & FAR.x<=1;
    B(1:sum(valid),c) = FAR.FB(valid);
    B(sum(valid):end,c) = B(sum(valid),c);
    PSS(1:sum(valid),c) = FAR.PSS(valid);
    PSS(sum(valid):end,c) = PSS(sum(valid),c);
    F(:,c) = k;
    c = c+1;
end

Then we set the colors for the colormap (which I partially took from you), and set the labels position:

colors = [0 0.2 0.4
    0 0.501 1;
    0 0.8 0.4;
    0.4 0.8 0;
    0.8 0.8 0;
    0.988 0.827 0.196;
    1 0.5019 0;
    0.8 0 0;
    0.6 0 0.2;
    0.4 0.1 0.5];

label_pos =[0.89 0.77
    1.01         0.74
    1.14         0.69
    1.37         0.64
    1.7          0.57
    2.03         0.41
    2.65         0.18
    2.925       -0.195
    2.75        -0.55];

And we use contourf to plot everything together, and set all kind of properties to make it look good:

[C,h] = contourf(B,PSS,F);
xlim([0 4])
ylim([-1 1])
colormap(colors)
caxis([0 1])
xlabel('Frequency Bias B')
ylabel('Pierce Skill Score PSS')
title('False Alarm Ratio')
ax = h.Parent;
ax.XTick = 0:4;
ax.YTick = -1:0.5:1;
ax.FontSize = 20;
for k = 1:numel(f)-2
    text(label_pos(k,1),label_pos(k,2),num2str(f(k+1)),...
        'FontSize',12+k)
end

And here is the result:

---

Getting the labels position:

If you wonder what is a fast way to obtain the variable label_pos, then here is how I made it...

You run the code above without the last for loop. Then you run the following code:

clabel(C,'manual')
f = gcf;
label_pos = zeros(numel(f.Children.Children)-1,2);
for k = 1:2:size(label_pos,1)
    label_pos(k,:) = f.Children.Children(k).Position(1:2);
end
label_pos(2:2:size(label_pos,1),:) = [];

After the first line the script will pause and you will see this message in the command window:

Carefully select contours for labeling. When done, press RETURN while the Graph window is the active window.

Click on the figure where you want to have a label, and press Enter.
That's it! Now the variable label_pos has the positions of the labels, just as I used it above.