You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
212 lines
6.9 KiB
Matlab
212 lines
6.9 KiB
Matlab
function [bss, rmse, dVc, dVm, ddV] = fsbeach_modelskill(measured, computed, initial)
|
|
%edited by Joshua Simmons 04/04/2016 to revise definitions of r2 etc
|
|
%XB_SKILL Computes a variety skill scores
|
|
%
|
|
% Computes a variety skill scores: R^2, Sci, Relative bias, Brier Skill
|
|
% Score. Special feature: within the XBeach testbed, the results are
|
|
% stored to be able to show the development of the different skill scores
|
|
% in time.
|
|
%
|
|
% Syntax:
|
|
% [r2 sci relbias bss] = xb_skill(measured, computed, varargin)
|
|
%
|
|
% Input:
|
|
% measured = Measured data where the first column contains independent
|
|
% values and the second column contains dependent values
|
|
% computed = Computed data where the first column contains independent
|
|
% values and the second column contains dependent values
|
|
% initial = Initial data where the first column contains independent
|
|
% values and the second column contains dependent values
|
|
% varargin = var: Name of the variable that is supplied
|
|
%
|
|
% Output:
|
|
% bss = Brier Skill Score
|
|
% rb20 = Regression through origin R2 value
|
|
%
|
|
% Example
|
|
% [r2 sci relbias bss] = xb_skill(measured, computed)
|
|
% [r2 sci relbias bss] = xb_skill(measured, computed, 'var', 'zb')
|
|
%
|
|
% See also xb_plot_skill
|
|
|
|
%% Copyright notice
|
|
% --------------------------------------------------------------------
|
|
% Copyright (C) 2011 Deltares
|
|
% Bas Hoonhout
|
|
%
|
|
% bas.hoonhout@deltares.nl
|
|
%
|
|
% Rotterdamseweg 185
|
|
% 2629HD Delft
|
|
%
|
|
% This library is free software: you can redistribute it and/or
|
|
% modify it under the terms of the GNU Lesser General Public
|
|
% License as published by the Free Software Foundation, either
|
|
% version 2.1 of the License, or (at your option) any later version.
|
|
%
|
|
% This library is distributed in the hope that it will be useful,
|
|
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
% Lesser General Public License for more details.
|
|
%
|
|
% You should have received a copy of the GNU Lesser General Public
|
|
% License along with this library. If not, see <http://www.gnu.org/licenses/>.
|
|
% --------------------------------------------------------------------
|
|
|
|
% This tool is part of <a href="http://OpenEarth.nl">OpenEarthTools</a>.
|
|
% OpenEarthTools is an online collaboration to share and manage data and
|
|
% programming tools in an open source, version controlled environment.
|
|
% Sign up to recieve regular updates of this function, and to contribute
|
|
% your own tools.
|
|
|
|
%% Version <http://svnbook.red-bean.com/en/1.5/svn.advanced.props.special.keywords.html>
|
|
% Created: 13 Apr 2011
|
|
% Created with Matlab version: 7.9.0.529 (R2009b)
|
|
|
|
% $Id: xb_skill.m 4725 2011-06-27 12:58:48Z hoonhout $
|
|
% $Date: 2011-06-27 14:58:48 +0200 (ma, 27 jun 2011) $
|
|
% $Author: hoonhout $
|
|
% $Revision: 4725 $
|
|
% $HeadURL: https://repos.deltares.nl/repos/OpenEarthTools/trunk/matlab/applications/xbeach/xb_analysis/xb_skill.m $
|
|
% $Keywords: $
|
|
|
|
%% read options
|
|
|
|
if ~exist('initial', 'var'); initial = []; end;
|
|
|
|
%% xbeach bugfix
|
|
% bug - sometime xbeach has suprious points on the end (a flat end to
|
|
% the profile) this is fake so remove it.
|
|
if length(computed)>=10
|
|
if nanmean(computed(1:10,2)) < nanmean(computed(end-9:end,2))
|
|
%xbeach direction
|
|
spr_len = length(diff(computed(:,2))) - find(diff(computed(:,2))~=0,1,'last');
|
|
computed(end-spr_len+1:end,2) = NaN;
|
|
else
|
|
spr_len = find(diff(computed(:,2))~=0,1,'first');
|
|
computed(1:spr_len-1,2) = NaN;
|
|
end
|
|
if spr_len > 2
|
|
disp(['Warning check computed points, ' num2str(spr_len) ' points found to be spurious'] )
|
|
end
|
|
end
|
|
|
|
%% remove nans
|
|
|
|
measured = measured(~any(isnan(measured), 2),:);
|
|
computed = computed(~any(isnan(computed), 2),:);
|
|
|
|
if ~isempty(initial)
|
|
initial = initial (~any(isnan(initial ), 2),:);
|
|
end
|
|
|
|
%% compute skills
|
|
|
|
[bss rmse dVc dVm ddV] = deal(nan);
|
|
|
|
if size(measured,1) > 1 && size(computed,1) > 1
|
|
|
|
x = unique([computed(:,1) ; measured(:,1)]);
|
|
x = x(~isnan(x));
|
|
|
|
zmt = interp1(measured(:,1), measured(:,2), x);
|
|
zct = interp1(computed(:,1), computed(:,2), x);
|
|
|
|
% determine active zone
|
|
if ~isempty(initial)
|
|
zit = interp1(initial(:,1), initial(:,2), x);
|
|
|
|
%average over 10 steps for small spikes
|
|
dzmi = smooth(abs(zmt-zit),10);
|
|
dzmi((isnan(zmt)|isnan(zit))) = 0;
|
|
dzcm = smooth(abs(zct-zmt),10);
|
|
dzcm((isnan(zct)|isnan(zmt))) = 0;
|
|
|
|
dz = max(dzmi,dzcm);
|
|
%has to be more change than a standard measurement error (10cm say)
|
|
zi1 = find(dz>0.1,1,'first');
|
|
zi2 = find(dz>0.1,1,'last');
|
|
|
|
x = x (zi1:zi2);
|
|
zmt = zmt(zi1:zi2);
|
|
zct = zct(zi1:zi2);
|
|
zit = zit(zi1:zi2);
|
|
end
|
|
|
|
%% clear all nan values in input (to avoid problems with interpolation)
|
|
%all equal at this point and assign values for BSS calc
|
|
x0 = x;
|
|
xc = x;
|
|
xm = x;
|
|
z0 = zit;
|
|
zc = zct;
|
|
zm = zmt;
|
|
|
|
x0(isnan(z0)) = [];
|
|
z0(isnan(z0)) = [];
|
|
xc(isnan(zc)) = [];
|
|
zc(isnan(zc)) = [];
|
|
xm(isnan(zm)) = [];
|
|
zm(isnan(zm)) = [];
|
|
|
|
%% determine new grid covered by all profiles
|
|
x_new = unique([x0' xc' xm']);
|
|
|
|
%% interpolate profiles onto new grid
|
|
[x0,id1] = unique(x0);
|
|
[xc,id2] = unique(xc);
|
|
[xm,id3] = unique(xm);
|
|
z0 = z0(id1);
|
|
zc = zc(id2);
|
|
zm = zm(id3);
|
|
|
|
z0_new = interp1(x0, z0, x_new);
|
|
zc_new = interp1(xc, zc, x_new);
|
|
zm_new = interp1(xm, zm, x_new);
|
|
|
|
% reduce vectors if any nans are left
|
|
id = ~(isnan(z0_new) | isnan(zc_new) | isnan(zm_new));
|
|
[x_new z0_new zc_new zm_new] = deal(x_new(id), z0_new(id), zc_new(id), zm_new(id));
|
|
|
|
%% Calculate volume change
|
|
try
|
|
Vc = volumefind(x_new,zc_new,0);
|
|
Vm = volumefind(x_new,zm_new,0);
|
|
V0 = volumefind(x_new,z0_new,0);
|
|
|
|
dVm = Vm-V0; %measured volume change
|
|
dVc = Vc-V0; %computed volume change
|
|
ddV = dVc-dVm; %predicted volume change - measured
|
|
nddV = ddV/abs(V0-Vm);
|
|
catch
|
|
dVm = NaN;
|
|
dVc = NaN;
|
|
ddV = NaN;
|
|
nddV = NaN;
|
|
end
|
|
%% derive weight
|
|
% weight is defined as half of diff left and half of diff right of each
|
|
% point; first and last point are considered to have only one side.
|
|
weight = diff([x_new(1) diff(x_new)/2 + x_new(1:end-1) x_new(end)]);
|
|
weight = weight / sum(weight);
|
|
|
|
%% calculate BSS only
|
|
% brier skill score
|
|
if ~isempty(initial)
|
|
mse_p = sum( (zm_new - zc_new).^2 .* weight );
|
|
mse_0 = sum( (zm_new - z0_new).^2 .* weight );
|
|
bss = 1. - (mse_p/mse_0);
|
|
else
|
|
bss = 1-(std(zc-zm))^2/var(zm);
|
|
end
|
|
|
|
%% Calculate others
|
|
|
|
rmse = sqrt(mean((zc_new-zm_new).^2));
|
|
rmsm = sqrt(mean(zm_new.^2));
|
|
sci = rmse/max(rmsm,abs(mean(zm_new)));
|
|
relbias = mean(zc_new-zm_new)/max(rmsm,abs(mean(zm_new)));
|
|
|
|
end
|
|
|