function varargout = debn_generate_wavedir_data(varargin)
%DEBN_GENERATE_WAVEDIR_DATA  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = debn_generate_wavedir_data(varargin)
%
%   Input: For <keyword,value> pairs call debn_generate_wavedir_data() without arguments.
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   debn_generate_wavedir_data
%
%   See also

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2012 Delft University of Technology
%       Kees den Heijer
%
%       C.denHeijer@TUDelft.nl
%
%       Faculty of Civil Engineering and Geosciences
%       P.O. Box 5048
%       2600 GA Delft
%       The Netherlands
%
%   This library is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 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 General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this library.  If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

% This tool is part of <a href="http://www.OpenEarth.eu">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: 19 Sep 2012
% Created with Matlab version: 7.12.0.635 (R2011a)

% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords: $

%%
OPT = struct(...
    'waterlevel_mainurl', 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/waterbase/sea_surface_height/',...
    'wavedir_mainurl', 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/waterbase/sea_surface_wave_from_direction/',...
    'destdir', fullfile(fileparts(which(mfilename)), '..', 'netica'));

OPT = setproperty(OPT, varargin);

%%
% % wave direction at IJmuiden munitiestortplaats
% ncfile = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/waterbase/sea_surface_wave_from_direction/id23-IJMDMNTSPS.nc';
% wavedir = nc_varget(ncfile, 'sea_surface_wave_from_direction');

%%
urls = opendap_catalog(OPT.waterlevel_mainurl);

for i = 1:length(urls)
    WL(i).url = urls{i};
    try
        WL(i).lon = nc_varget(urls{i}, 'lon');
        WL(i).lat = nc_varget(urls{i}, 'lat');
        [WL(i).x,WL(i).y] = convertCoordinates(WL(i).lon,WL(i).lat,'CS1.code',4326,'CS2.code',28992);
%         WL(i).x = nc_varget(urls{i}, 'x');
%         WL(i).y = nc_varget(urls{i}, 'y');
    end
end

id = ~cellfun(@isempty,{WL.y}) & ~cellfun(@isempty,{WL.x});
WL = WL(id);

%%
urls = opendap_catalog(OPT.wavedir_mainurl);

for i = 1:length(urls)
    A(i).url = urls{i};
    try
        A(i).lon = nc_varget(urls{i}, 'lon');
        A(i).lat = nc_varget(urls{i}, 'lat');
        [A(i).x,A(i).y] = convertCoordinates(A(i).lon,A(i).lat,'CS1.code',4326,'CS2.code',28992);
%         A(i).x = nc_varget(urls{i}, 'x');
%         A(i).y = nc_varget(urls{i}, 'y');
    end
end

% % distance from individual wave angle locations to all water level stations
% D = cellfun(@(x,y) sqrt((x-[WL.x]).^2+(y-[WL.y]).^2), {A.x}, {A.y},...
%     'uniformoutput', false);

ixA = ismember([A.lat], [WL.lat]) & ismember([A.lon], [WL.lon]);
A = A(ixA);
ixWL = ismember([WL.lat], [A.lat]) & ismember([WL.lon], [A.lon]);
WL = WL(ixWL);

for i = 1:length(WL)
    A(i) = A(WL(i).x == [A.x]);
end

% % identify nearest water level stations
% [v ixWL] = cellfun(@min, D);
Awl = struct('WLurl', {WL.url});
A = mergestructs(A, Awl);

[d, h, t] = cellfun(@commontime, {A.url}, {A.WLurl}, 'uniformoutput', false);

[~, ix] = cellfun(@(x) sort(x, 'descend'), h, 'uniformoutput', false);

[ds, hs] = cellfun(@(x,y,ix) deal(x(ix), y(ix)), d, h, ix, 'uniformoutput', false);

wavedir = cell2mat(cellfun(@(x) x(1:round(length(x)/500))', ds, 'uniformoutput', false));

%%
netica_write_case_file(fullfile(OPT.destdir, 'wavedir.cas'),...
    'wavedir', mod(wavedir, 360));

%%
ncfile = 'http://opendap.deltares.nl/thredds/dodsC/opendap/deltares/landboundaries/northsea.nc';
x = nc_varget(ncfile, 'x');
y = nc_varget(ncfile, 'y');

[~, ix] = unique([A.y]); % unique automatically sorts ascending
ix = fliplr(ix); % flip to sort descending (north to south)
stations = [{A(ix).x}; {A(ix).y}; repmat({'o'}, size(ix))];
stationnames = cellfun(@(ncfile) nc_varget(ncfile, 'station_name'), {A(ix).url}, 'uniformoutput', false);

figure;
ph = plot(x,y,...
    stations{:});
set(ph(1), 'HandleVisibility', 'off');
set(ph(2:end), {'displayname'}, stationnames(:))

axis image
xlabel('x coordinate [m]')
ylabel('y coordinate [m]')

lh = legend(gca,...
    'location', 'northwest');

figureStyle(@styleThesis,...
    {'figure' 'Position' [1 1 15 10]});
tightenfigure;
[fpath fname] = fileparts(which(mfilename));
% fname = sprintf('%s_%s', fname, strrep(OPT.suffix, ' ', '_'));
figurepath = fullfile(fpath, '..', 'figures');
set(gcf, 'tag', '$HeadURL$ $Revision$')
print(gcf, '-depsc2', fullfile(figurepath, [fname '.eps']))
if ispc
    saveas(gcf, fullfile(figurepath, [fname '.emf']))
end
saveas(gcf, fullfile(figurepath, [fname '.fig']))

function [d, h, t] = commontime(url1, url2)

t1 = nc_varget(url1, 'time');
t2 = nc_varget(url2, 'time');
d1 = nc_varget(url1, 'sea_surface_wave_from_direction');
h2 = nc_varget(url2, 'sea_surface_height');

[t, ix1, ix2] = intersect(t1, t2);
d = d1(ix1);
h = h2(ix2);
