data_dir = 'C:\OpenEarthRawData\rijkswaterstaat\strandlijnen\raw\data\';
files = dir([data_dir '*98']);
%% read XLS files
data = [];
for ii = 1:length(files);
    disp(['processing ' files(ii).name]);
    fid = fopen([data_dir files(ii).name]);
    C = textscan(fid,'%f%f%f%f%f','CollectOutput',1);
    fclose(fid);
    switch files(ii).name
        case 'Str_schi_98', kustvak =  2;
        case 'Str_amel_98', kustvak =  3;
        case 'Str_ters_98', kustvak =  4;
        case 'Str_vlie_98', kustvak =  5;
        case 'Str_texe_98', kustvak =  6;
        case 'Str_nho_98' , kustvak =  7;
        case 'Str_rijn_98', kustvak =  8;
        case 'Str_delf_98', kustvak =  9;

        case 'Str_voor_98', kustvak = 11;
        case 'Str_goer_98', kustvak = 12;
        case 'Str_scho_98', kustvak = 13;

        case 'Str_walc_98', kustvak = 16;
        case 'Str_zvl_98' , kustvak = 17;
    end
    data = [data;[repmat(kustvak,size(C{1,1}(:,1))) C{1,1}]];
end
%% manual adjustments 
% remove a non existing transect:
% transect 5600 in rijnland, this transect is in tghe IJmuiden harbour
ind = find(data(:,1) == 8 & data(:,2)==56);
data(ind,:) = [];
ind = find(data(:,1) == 13 & data(:,2)==9.32);
data(ind,2) = 9.39;

% replace fill value with nan
data(data==99999)=NaN;

%% some data processing
R.kustvak  = data(:,1);
R.transect = data(:,2)*1000;
R.year     = data(:,3);
R.DF       = data(:,4);
R.MHW      = data(:,5);
R.MLW      = data(:,6);
R.trID     = R.kustvak*1e6+R.transect;

%% reshape the data
% pre allocating
D.year = unique(R.year);
D.trID = unique(R.trID);

D.kustvak  = nan(length(D.trID),length(D.year));
D.transect = nan(length(D.trID),length(D.year));
D.DF       = nan(length(D.trID),length(D.year));
D.MHW      = nan(length(D.trID),length(D.year));
D.MLW      = nan(length(D.trID),length(D.year));


% loop through the data
for jj = 1:length(D.year);
    disp(['reshaping ' num2str(D.year(jj))]);
    for ii = 1:length(D.trID)
        ind = find(R.trID==D.trID(ii)&R.year==D.year(jj));
        if ~isempty(ind)
            D.kustvak(ii,jj)  = R.kustvak(ind);
            D.transect(ii,jj) = R.transect(ind);
            D.DF(ii,jj)       = R.DF(ind);
            D.MHW(ii,jj)      = R.MHW(ind);
            D.MLW(ii,jj)      = R.MLW(ind);
        end
    end
end

%% convert coordinates
EPSG = load('EPSGnew');
disp('converting coordinates...');
[xRD,yRD] = xRSP2xyRD(D.DF(:),D.kustvak(:),D.transect(:));
[lon,lat] = convertCoordinates(xRD,yRD,EPSG,'CS1.name','Amersfoort / RD New','CS2.name','WGS 84','CS2.type','geo');
D.xDF   = reshape(xRD,size(D.DF));
D.yDF   = reshape(yRD,size(D.DF));
D.latDF = reshape(lat,size(D.DF));
D.lonDF = reshape(lon,size(D.DF));

disp('converting more coordinates...');
[xRD,yRD] = xRSP2xyRD(D.MHW(:),D.kustvak(:),D.transect(:));
[lon,lat] = convertCoordinates(xRD,yRD,EPSG,'CS1.name','Amersfoort / RD New','CS2.name','WGS 84','CS2.type','geo');
D.xMHW   = reshape(xRD,size(D.DF));
D.yMHW   = reshape(yRD,size(D.DF));
D.latMHW = reshape(lat,size(D.DF));
D.lonMHW = reshape(lon,size(D.DF));

disp('almost done...');
[xRD,yRD] = xRSP2xyRD(D.MLW(:),D.kustvak(:),D.transect(:));
[lon,lat] = convertCoordinates(xRD,yRD,EPSG,'CS1.name','Amersfoort / RD New','CS2.name','WGS 84','CS2.type','geo');
D.xMLW   = reshape(xRD,size(D.DF));
D.yMLW   = reshape(yRD,size(D.DF));
D.latMLW = reshape(lat,size(D.DF));
D.lonMLW = reshape(lon,size(D.DF));

%% Save D
save('StrandLijnData.mat','D');
disp('done');