% compareMW
%
% compare mean high water and mean low water of the old database and the
% new one, based on the intersection of the Jarkus profiles (until 2012)..
% Therefore the comparison is performed just on the period 1965 - 1998.
%
% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords$ 

clear 
close all

% old database: standlijnen.nc for info meanhigwater, meanlowwater, dunefoot 1843-1998
OPT.ncpath3   = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/strandlijnen/';
OPT.Filename3 = 'strandlijnen.nc';
% read from 1843-1998 strandlijnen.nc file
time_STRAND = str2num(datestr(nc_varget([OPT.ncpath3 OPT.Filename3],'time') + datenum('1970','yyyy'),'yyyy')); % year
MHW_STRAND  = nc_varget([OPT.ncpath3 OPT.Filename3],'MHW');
MLW_STRAND  = nc_varget([OPT.ncpath3 OPT.Filename3],'MLW');
DF_STRAND   = nc_varget([OPT.ncpath3 OPT.Filename3],'DF');
id_STRAND   = nc_varget([OPT.ncpath3 OPT.Filename3],'id');

% new database: jarkus database - dataset 1965
OPT.ncpath2      = 'c:\\Users\santinel\Documents\openearthrawdata\rijkswaterstaat\jarkus\scripts\';
OPT.Filename2    = 'transect.nc';
% Read data
id       = nc_varget([OPT.ncpath2 OPT.Filename2],'id');
time     = str2num(datestr(nc_varget([OPT.ncpath2 OPT.Filename2],'time') + datenum('1970','yyyy'),'yyyy')); % year
areacode = nc_varget([OPT.ncpath2 OPT.Filename2],'areacode');
cross    = nc_varget([OPT.ncpath2 OPT.Filename2],'cross_shore');
mean_high_water = nc_varget([OPT.ncpath2 OPT.Filename2],'mean_high_water');
mean_low_water  = nc_varget([OPT.ncpath2 OPT.Filename2],'mean_low_water');

[idIDa idIDb]       = ismember(id,id_STRAND); % 324 common ones of 334 from strandlijnen
[timeIDa timeIDb]   = ismember(time,time_STRAND);

% read from .mat file created with MW_Jarkus.m
load('MHW.mat');
load('MLW.mat');

% counters
[kv, m, dummy] = unique(areacode);
n          = vertcat(1, m(1:end-1)+1);

% make dir
if ~exist('Comparison'); mkdir('Comparison'); end

for ii = 1:length(id) 
    %num_tr  = m(i)-n(i)+1;
    if idIDa(ii)
        altitude = nc_varget([OPT.ncpath2 OPT.Filename2],'altitude',[0 ii-1 0],[-1 1 -1]);
        disp(['comparison: transect ' num2str(ii) ' of ' num2str(length(id))]);
        
        for k = 1:sum(timeIDa) % loop for years in common
            if ~isempty(altitude(k,~isnan(altitude(k,:)))) % altitude measurement exists that year
                figure('Visible', 'Off')
                set(gcf,'Position',[5 35 1040 760]);
                
                % plot profile
                plot((cross(~isnan(altitude(k,:))))',altitude(k,~isnan(altitude(k,:))), 'Color', [1 .5 0], 'LineWidth',2);
                hold on; grid on;
                % plot the values of MHW MLW from 1990
                L1 = plot(get(gca,'XLim'), [mean_high_water(ii), mean_high_water(ii)], 'Color', [0 .5 0], 'LineWidth',2);
                plot(get(gca,'XLim'), [mean_low_water(ii), mean_low_water(ii)], 'Color', [0 .5 0], 'LineWidth',2);
                % plot MHW, MLW strandlijnen (RED); 
                L2 = plot([MHW_STRAND(idIDb(ii),timeIDb(k)), MHW_STRAND(idIDb(ii),timeIDb(k))], get(gca,'YLim'), 'Color', [1 0 0], 'LineWidth',1);
                plot([MLW_STRAND(idIDb(ii),timeIDb(k)), MLW_STRAND(idIDb(ii),timeIDb(k))], get(gca,'YLim'), 'Color', [1 0 0], 'LineWidth',1);
                % plot MHW, MLW from transect (BLUE);
                L3 = plot([MHW_cross(ii,k), MHW_cross(ii,k)], get(gca,'YLim'), 'Color', [0 0 1], 'LineWidth',1);
                plot([MLW_cross(ii,k), MLW_cross(ii,k)], get(gca,'yLim'), 'Color', [0 0 1], 'LineWidth',1);
                
                title(['Comparison: transect ' num2str(id(ii))  ', year ' num2str(time(k))], 'FontSize', 12)
                legend([L1, L2, L3],{'MeanWater90 ','MeanWaterold','MeanWaternew'}, 'location', 'NorthEast');
                print('-dpng','-r300',['Comparison' filesep 'transect' num2str(id(ii)) '_year' num2str(time(k))]);
                close
            
            end
        end
    end

end
