function [seriesTime, seriesTimeHeightData, seriesWaveHeight] = rws_waterbase_waveclimate(varargin)
%RWS_WATERBASE_WAVECLIMATE Creates wave climate from Rijkswaterstaat waterbase data.
%
% This function takes the netcdf version from the data from the waterbase website from
% Rijkswaterstaat, and creates a wave climate for certain locations and
% time period, weighted by wave height. If no timeframe is provided, a GUI will pop up.
%
% The required data can be obtained from the Deltares
% OpenDaP server (default argument). Note however that these data are
% outdated. To obtain the latest data, create a local netcdf copy using
% rws_waterbase_all.
%
% Syntax:
% varargout = rws_waterbase_waveclimate(varargin)
%
% Input:
% varargin =
% 'baseurl' = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/waterbase/';
% 'time' = []; % Start and end date of wave input data. Format: [yyyymmdd yyyymmdd]
% 'outputrawdata' = 1; % output raw data to file?
% 'orawfilename' = ''; % filename for output raw data, leave empty if multiple stations
% 'classificationpars' = []; % Vector with classif. params. Format: [hs0,dhs,nhs,dir0,ddir,ndir]
% 'oclimfilename' = ''; % filename for output climate data, leave empty if multiple stations
%
%
% Output:
% varargout =
% none
%
% Example
% rws_waterbase_waveclimate({'outputrawdata',0,'classificationpars',[0 0.5 5 5 10 36]})
%
% See also
% rws_waterbase_get_url, rws_waterbase_get_locations,
% rws_waterbase_all
%% Copyright notice
% --------------------------------------------------------------------
% Copyright (C) 2013 UNESCO-IHE
% Johan Reyns, Dano Roelvink
%
% {j.reyns,d.roelvink}@unesco-ihe.org
%
% P.O. Box 3015
% 2601 DA Delft
%
% 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 .
% --------------------------------------------------------------------
% This tool is part of OpenEarthTools.
% 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
% Created: 15 Apr 2013
% Created with Matlab version: 7.8.0.347 (R2009a)
% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords: $
%% Init
OPT.strmatch = 'exact';
OPT.version = 2; % 1 is before summer 2009, 2 is after mid dec 2009
OPT.baseurl = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/waterbase/';
OPT.time = []; % dummy
OPT.outputrawdata = 1; % output raw data to file?
OPT.orawfilename = ''; % filename for output raw data, leave empty if multiple stations
OPT.classificationpars = []; % hs0,dhs,nhs,dir0,ddir,ndir
OPT.oclimfilename = ''; % filename for output climate data, leave empty if multiple stations
OPT = setproperty(OPT,varargin{:});
%% Get pre-defined location names
LOC = [];
while isempty(LOC)
LOC = rws_waterbase_get_locations(22,'22%7CSignificante+golfhoogte+uit+energiespectrum+van+30-500+mhz+in+cm+in+oppervlaktewater');
if isempty(LOC)
fprintf(2,[OPT.baseurl,' failed. Trying again to get list of: Stations']);
pause(5)
end
end
[indLoc, ok] = listdlg('ListString', LOC.FullName, ...
'SelectionMode', 'multiple', ...
'InitialValue' , [1:length(LOC.FullName)], ...
'PromptString' , 'Select the locations', ....
'Name' , 'Selection of locations');
if (ok == 0); return; end
if length(indLoc)>1
disp(['message: rws_waterbase_waveclimate: Location ',num2str(length(indLoc),'%0.3d'),'x: #',num2str(indLoc(:)','%0.3d,')])
else
disp(['message: rws_waterbase_waveclimate: Location # ',num2str(indLoc,'%0.3d'),': ',LOC.ID{indLoc},' "',LOC.FullName{indLoc},'"'])
end
%% Times
if ~isempty(OPT.time)
indDate = OPT.time;
startdate = [datestr(indDate(1),'yyyymmddHHMM')]; %,'01010000'];
enddate = [datestr(indDate(2),'yyyymmddHHMM')]; %,'12312359'];
ok = 1;
else
ListYear = '1961';
for iYear = 1962:str2num(datestr(date,'yyyy'))
ListYear = strvcat(ListYear, sprintf('%d', iYear));
end
ListYear = cellstr(ListYear);
[indDate, ok] = listdlg('ListString', ListYear, ...
'SelectionMode', 'multiple', ...
'InitialValue' , [length(ListYear)], ...
'PromptString' , 'Select the year', ....
'Name' , 'Selection of year');
if (ok == 0)
return;
end
startdate = [ListYear{min(indDate)} '01010000'];
enddate = [ListYear{max(indDate)} '12312359'];
end
disp(['message: rws_waterbase_waveclimate: startdate ',startdate]);
disp(['message: rws_waterbase_waveclimate: enddate ',enddate]);
h = waitbar(0,'Downloading data...');
for iLoc=1:length(indLoc)
urlNameWaveHeight = [OPT.baseurl,'17_Significante_golfhoogte_uit_energiespectrum_van_30-500_mhz_in_cm_in_oppervlaktewater/nc/','id22-',LOC.ID{indLoc(iLoc)},'.nc'];
urlNameWaveDirection = [OPT.baseurl,'08_Gemiddelde_richting_uit_golfrichtingsspectrum_van_30-500_mhz_in_graad_t.o.v._kaart_noorden_in_oppervlaktewater/nc/' ,'id23-',LOC.ID{indLoc(iLoc)},'.nc'];
urlNameWavePeriod = [OPT.baseurl,'07_Gem._golfperiode_uit_spectrale_momenten_m0+m2_van_30-500_mhz_in_s_in_oppervlaktewater/nc/' ,'id24-',LOC.ID{indLoc(iLoc)},'.nc'];
%% Perform operation on three separate time series, as they are not
%% necessarily synchronous
seriesTimeHeight = nc_varget(urlNameWaveHeight,'time');
seriesTimePeriod = nc_varget(urlNameWavePeriod,'time');
seriesTimeDirection = nc_varget(urlNameWaveDirection,'time');
seriesTimeHeightStart = find(datenum('010119700100','ddmmyyyyHHMM')+seriesTimeHeight >=datenum(startdate,'yyyymmddHHMM'),1,'first');
seriesTimeHeightEnd = find(datenum('010119700100','ddmmyyyyHHMM')+seriesTimeHeight <=datenum(enddate ,'yyyymmddHHMM'),1,'last');
seriesTimePeriodStart = find(datenum('010119700100','ddmmyyyyHHMM')+seriesTimePeriod >=datenum(startdate,'yyyymmddHHMM'),1,'first');
seriesTimePeriodEnd = find(datenum('010119700100','ddmmyyyyHHMM')+seriesTimePeriod <=datenum(enddate ,'yyyymmddHHMM'),1,'last');
seriesTimeDirStart = find(datenum('010119700100','ddmmyyyyHHMM')+seriesTimeDirection>=datenum(startdate,'yyyymmddHHMM'),1,'first');
seriesTimeDirEnd = find(datenum('010119700100','ddmmyyyyHHMM')+seriesTimeDirection<=datenum(enddate ,'yyyymmddHHMM'),1,'last');
if numel(seriesTimeHeightStart)==0 || numel(seriesTimeHeightEnd)==0
disp(['Time limits outside period of available data.. Data limits are ', datestr(datenum('01011970','ddmmyyyy')+seriesTimeHeight(1)), ' and ',datestr(datenum('01011970','ddmmyyyy')+seriesTimeHeight(end)),'. Aborting.']);
close(h);
return;
end
seriesTimeHeightData = datenum('010119700100','ddmmyyyyHHMM')+ seriesTimeHeight(seriesTimeHeightStart:seriesTimeHeightEnd)';
seriesTimePeriodData = datenum('010119700100','ddmmyyyyHHMM')+ seriesTimePeriod(seriesTimePeriodStart:seriesTimePeriodEnd)';
seriesTimeDirectionData = datenum('010119700100','ddmmyyyyHHMM')+ seriesTimeDirection(seriesTimeDirStart:seriesTimeDirEnd )';
seriesWaveHeight = nc_varget(urlNameWaveHeight ,'sea_surface_wave_significant_height',[0 seriesTimeHeightStart-1],[1 seriesTimeHeightEnd-seriesTimeHeightStart+1]);
seriesWaveDir = nc_varget(urlNameWaveDirection,'sea_surface_wave_from_direction' ,[0 seriesTimeDirStart-1],[1 seriesTimeDirEnd-seriesTimeDirStart+1]);
seriesWavePeriod = nc_varget(urlNameWavePeriod ,'sea_surface_wind_wave_tm02' ,[0 seriesTimePeriodStart-1],[1 seriesTimePeriodEnd-seriesTimePeriodStart+1]); % Tm02
resolution = seriesTimeHeightData(3)-seriesTimeHeightData(2);
resolution = (ceil(resolution*24*60))/24/60;
%% Plot raw data
figure(iLoc*2-1);
subplot(3,1,1);
plot(seriesTimeHeightData,seriesWaveHeight,'-r');
ylabel('Wave height [m]')
datetick('x')
subplot(3,1,2);
plot(seriesTimePeriodData,seriesWavePeriod,'-g');
ylabel('Wave period [s]')
datetick('x')
subplot(3,1,3);
ylim([0 360]);
plot(seriesTimeDirectionData,seriesWaveDir,'-b');
datetick('x');
ylabel('Wave direction [deg from N]')
xlabel('time')
%% Generate synchronous time series
seriesTime = (floor(datenum(startdate,'yyyymmddHHMM')):resolution:ceil(datenum(enddate,'yyyymmddHHMM'))-1/24)';
seriesWaveHeightC = NaN(length(seriesTime),1);
seriesWavePeriodC = NaN(length(seriesTime),1);
seriesWaveDirC = NaN(length(seriesTime),1);
%% Fill series with data
[dummy,iWHa,iWHb] = intersect(round(seriesTime*10000)/10000,round(seriesTimeHeightData(:)*10000)/10000); % to counter rounding error related to ismember implementation
[dummy,iWPa,iWPb] = intersect(round(seriesTime*10000)/10000,round(seriesTimePeriodData(:)*10000)/10000);
[dummy,iWDa,iWDb] = intersect(round(seriesTime*10000)/10000,round(seriesTimeDirectionData(:)*10000)/10000);
seriesWaveHeightC(iWHa) = seriesWaveHeight(iWHb);
seriesWavePeriodC(iWPa) = seriesWavePeriod(iWPb);
seriesWaveDirC(iWDa) = seriesWaveDir (iWDb);
clear dummy;
%% Write raw data to file
if OPT.outputrawdata
if isempty(OPT.orawfilename)
OPT.orawfilename = ['wavedata_',LOC.ID{indLoc(iLoc)},'_',startdate,'_',enddate,'.txt'];
end
clear ii;
ofid = fopen(OPT.orawfilename,'w');
fprintf(ofid,'%s\t%s\t%s\t%s\n','Time [yyyymmddHHMM]','Hs [m]','Tp [s]','Dir [deg]');
for ii = 1:length(seriesTime)
fprintf(ofid,'%s\t%.2f\t%.2f\t%.2f\n',datestr(seriesTime(ii),'yyyymmddHHMM'),seriesWaveHeightC(ii),seriesWavePeriodC(ii),seriesWaveDirC(ii));
end
fclose(ofid); clear ofid;
OPT.orawfilename='';
end
%% Convert timeseries of wave height, direction and period into
% - 2D table of percentage occurrence of wave height classes perc
% - 2D table of weighted mean waveheight per class meanhs
% - 2D table of weighted mean direction per class meandir
% - 2D table of weighted mean wave period per class meantp
% wave heights are weighted with power 2.5
hs = seriesWaveHeightC(:);
dir = seriesWaveDirC(:);
tp = 1.2859*seriesWavePeriodC(:); % assuming JONSWAP, gamma=3.3
hs1 = hs(find(~isnan(hs)&~isnan(dir)&~isnan(tp)));
tp1 = tp(find(~isnan(hs)&~isnan(dir)&~isnan(tp)));
dir = dir(find(~isnan(hs)&~isnan(dir)&~isnan(tp)));
hs = hs1;
tp = tp1;
power=2.5;
if isempty(OPT.classificationpars)
hs0 = 0;
dhs = 0.5;
nhs = 5;
dir0 = 5;
ddir = 10;
ndir = 36;
else
if length(OPT.classificationpars)~=6
disp('Wrong number of arguments for wave climate classification. Aborting..');
clear all; fclose all;
return;
end
hs0 = OPT.classificationpars(1);
dhs = OPT.classificationpars(2);
nhs = OPT.classificationpars(3);
dir0 = OPT.classificationpars(4);
ddir = OPT.classificationpars(5);
ndir = OPT.classificationpars(6);
end
ind_hs =max(min(ceil((hs-hs0)/dhs),nhs),1);
ind_dir=max(min(ceil((dir-dir0)/ddir),ndir),1);
n = zeros(nhs,ndir);
s_hs = zeros(nhs,ndir);
s_dir = zeros(nhs,ndir);
s_tp = zeros(nhs,ndir);
for i=1:length(hs);
ih=ind_hs(i);
id=ind_dir(i);
n (ih,id) = n(ih,id)+1;
s_hs (ih,id) = s_hs (ih,id)+hs(i).^power;
s_dir(ih,id) = s_dir(ih,id)+dir(i);
s_tp (ih,id) = s_tp (ih,id)+tp(i);
end;
n(n==0)= NaN;
perc = n/length(hs)*100;
meanhs = (s_hs./n).^(1/power);
meandir = s_dir./n;
meantp = s_tp./n;
%% Write output p,Hs,Tp,Dir
if isempty(OPT.oclimfilename)
OPT.oclimfilename = ['waveclimate_',LOC.ID{indLoc(iLoc)},'_',startdate,'_',enddate,'.txt'];
end
ofid = fopen(OPT.oclimfilename,'w');
clear ii;
for ii = 1:size(perc,2)
for ij = 1:size(perc,1)
fprintf(ofid,'%.2f\t%.2f\t%.2f\t%.0f\n',perc(ij,ii),meanhs(ij,ii),meantp(ij,ii),meandir(ij,ii));
end
end
fclose(ofid); clear ofid;
OPT.oclimfilename='';
%% Plot output in spiderplot
figure(iLoc*2);
xhs=[hs0:dhs:nhs*dhs];
xdir=[dir0:ddir:ndir*ddir];
perc_ext=zeros([size(perc,1)+1,size(perc,2)+1]);
perc_ext(1:end-1,1:end-1)=perc;
meanhs_ext=zeros([size(meanhs,1)+1,size(meanhs,2)+1]);
meanhs_ext(1:end-1,1:end-1)=meanhs;
meantp_ext=zeros([size(meantp,1)+1,size(meantp,2)+1]);
meantp_ext(1:end-1,1:end-1)=meantp;
for i=1:ndir+1;
theta=(i-1)*ddir*pi/180;
for j=1:length(xhs);
xpol(j,i)=xhs(j)*sin(theta);
ypol(j,i)=xhs(j)*cos(theta);
end;
end
ticks=[-xhs(end) 0 xhs(end)];
ticklabels=[xhs(end) 0 xhs(end)];
subplot(131);
pcolor(xpol,ypol,perc_ext);
axis equal tight
colorbar;
set(gca,'fontweight','bold');
title('percentage occurrence');
set(gca,'ytick',ticks,'yticklabel',ticklabels,'xticklabel',[]);
ylabel('Hs (m)')
subplot(132);
pcolor(xpol,ypol,meanhs_ext);
axis equal tight
colorbar;
set(gca,'fontweight','bold');
title('mean Hs [m]');
set(gca,'xticklabel',[],'yticklabel',[]);
subplot(133);
pcolor(xpol,ypol,meantp_ext);
axis equal tight
colorbar;
set(gca,'fontweight','bold');
title('mean T_p [s]');
set(gca,'xtick',ticks,'xticklabel',ticklabels,'ytick',ticks,'yticklabel',ticklabels);
xlabel('Hs [m]');
ylabel('Hs [m]');
waitbar(iLoc/length(indLoc),h);
end
close(h);
varargout = {};
end