%% set oetsettings
close all; clear all;



%% initialize superTrans
addpath(genpath('~/Documents/checkouts/McTools/mc_programs/DelftAlmighty/SuperTrans'));
load('CoordinateSystems.mat')
load('Operations.mat');

%% settings
data_source = [1 0];     % data_source(1) = griddata and data_source(2) = transect data
d3d = 0;                 % d3d = 1, means Delft3D grid and bathymetry files are generated
XBeach = 1;              % XBeach = 1, means XBeach grid and bathymetry files are generated


%% set D3D toolboxpath
if d3d
    addpath F:\d3d\toolbox\
    wlsettings('F:\d3d\toolbox\');
end

%% griddata

if data_source(1)==1
    % TO DO: make generic function for easy use in future...
    % [X,Y,Z] = GetGridData(P)
    urls = {'http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/vaklodingen/KB127_1312.nc',
            'http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/vaklodingen/KB128_1312.nc',
            'http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/vaklodingen/KB129_1312.nc',
            'http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/vaklodingen/KB130_1312.nc'};
    urls={'/Users/fedorbaart/Documents/checkouts/OpenEarthRawData/trunk/rijkswaterstaat/vaklodingen/scripts/KB127_1312.nc',
          '/Users/fedorbaart/Documents/checkouts/OpenEarthRawData/trunk/rijkswaterstaat/vaklodingen/scripts/KB128_1312.nc',
          '/Users/fedorbaart/Documents/checkouts/OpenEarthRawData/trunk/rijkswaterstaat/vaklodingen/scripts/KB129_1312.nc',
          '/Users/fedorbaart/Documents/checkouts/OpenEarthRawData/trunk/rijkswaterstaat/vaklodingen/scripts/KB130_1312.nc'};
    urls={'/Users/fedorbaart/Documents/checkouts/OpenEarthRawData/trunk/rijkswaterstaat/vaklodingen/scripts/jarkusKB127_1312.nc',
          '/Users/fedorbaart/Documents/checkouts/OpenEarthRawData/trunk/rijkswaterstaat/vaklodingen/scripts/jarkusKB128_1312.nc',
          '/Users/fedorbaart/Documents/checkouts/OpenEarthRawData/trunk/rijkswaterstaat/vaklodingen/scripts/jarkusKB129_1312.nc',
          '/Users/fedorbaart/Documents/checkouts/OpenEarthRawData/trunk/rijkswaterstaat/vaklodingen/scripts/jarkusKB130_1312.nc'
          '/Users/fedorbaart/Documents/checkouts/OpenEarthRawData/trunk/rijkswaterstaat/vaklodingen/scripts/vaklodingenKB127_1312.nc',
          '/Users/fedorbaart/Documents/checkouts/OpenEarthRawData/trunk/rijkswaterstaat/vaklodingen/scripts/vaklodingenKB128_1312.nc',
          '/Users/fedorbaart/Documents/checkouts/OpenEarthRawData/trunk/rijkswaterstaat/vaklodingen/scripts/vaklodingenKB129_1312.nc',
          '/Users/fedorbaart/Documents/checkouts/OpenEarthRawData/trunk/rijkswaterstaat/vaklodingen/scripts/vaklodingenKB130_1312.nc'};
    x = []; y = []; Z = [];    
    Z_hist = [];
    enddate = 0;
    for i = 1:length(urls)
        % we have:
        % - old values (the ones we collected so far)
        % - extra values (the ones which we read in now)
        % - new values (the combination of old values and extra values)
        
        % store the old values
        old_x = x;
        old_y = y;
        old_Z = Z;
        old_Z_hist = Z_hist;
        
        url = urls{i};
        % info = nc_info(url);
        time = nc_varget(url,'time'); % Fedor checks whether database
        % date corresponds to matlab date
        
        % read extra values
        extra_x = nc_varget(url,'x');
        extra_y = nc_varget(url,'y');
        extra_Z = nc_varget(url,'z',[0,0,0],[1,-1,-1]); %read first time
        extra_Z_hist = zeros(size(extra_Z))*nan; % store the history
        for j=1:length(time)
            if enddate && time(j)+datenum('1-1-1') > enddate
                disp(time(j))
                break
            end
            % we create a temporary z variable because we want to store the
            % z's on top of each other
            tmp_z = nc_varget(url,'z',[j-1,0,0],[1,-1,-1]);
            disp(['Adding', datestr(time(j)+datenum('1-1-1')) 'for' url])
            extra_Z(~(isnan(tmp_z))) = tmp_z(~isnan(tmp_z)); % replace all values by newer values if not nan
            extra_Z_hist(~isnan(tmp_z)) = time(j); % store the date of each Z point
        end
        
        %% combine the old and extra values into a new value
        new_x = [old_x; extra_x]; % Fedor is going to add information to database to automatically merge the datasets
        new_x = sort(unique(new_x));
        new_y = [old_y; extra_y];
        new_y = sort(unique(new_y));
        
        % Putting the extra_Z values in the correct place is a bit
        % difficult
        % We also keep a history here so we know how old the data is. 
        % create a new Z of the correct size
        new_Z = zeros(length(new_y),length(new_x))*nan;
        new_Z_hist = zeros(size(new_Z))*nan;
        % The grid might have changed so we need to lookup where 
        % to put the old data
        [c,iold,inew] = intersect(old_y, new_y);
        [c,jold,jnew] = intersect(old_x, new_x);
        % put the data in these locations and store the history on the same
        % place
        new_Z(inew, jnew) = old_Z(iold,jold);
        new_Z_hist(inew, jnew) = old_Z_hist(iold,jold);
        % now we gonna add the extra data, again, we might have a different
        % grid. So lookup where the values should go.
        [c,iextra,inew] = intersect(extra_y, new_y);
        [c,jextra,jnew] = intersect(extra_x, new_x);
        % We know where to put the extra value but we only want to store
        % values if they are not nan.
        % we first extract the place where we want to store the data
        % this is required because we can't say:
        % new_Z(inew,jnew)(~isnan(extra_Z(iextra,jextra)))=
        % extra_Z(iextra,jextra)(~isnan(extra_Z(iextra,jextra))
        % now we have to put it in 8 lines
        
%         Inew = false(size(new_Z));
%         [Inew(inew,jnew), Iextra] = deal(~isnan(extra_Z(iextra,jextra)));
%         new_Z(Inew) = extra_Z(Iextra);
%         new_Z_hist(Inew) = extra_Z_hist(Iextra);

        slice_Z = new_Z(inew,jnew);
        slice_Z_hist = new_Z_hist(inew,jnew);
        extra_slice_Z = extra_Z(iextra,jextra);
        extra_slice_Z_hist = extra_Z_hist(iextra, jextra);
        slice_Z(~isnan(extra_slice_Z)) = extra_slice_Z(~isnan(extra_slice_Z));
        slice_Z_hist(~isnan(extra_slice_Z)) = extra_slice_Z_hist(~isnan(extra_slice_Z));
        new_Z(inew,jnew) = slice_Z;
        new_Z_hist(inew,jnew) = slice_Z_hist;
        
        % store the new arrays
        x = new_x;
        y = new_y;
        Z = new_Z;
        Z_hist = new_Z_hist;
       
    end
    [XV,YV] = meshgrid(x,y);
    ZV = Z;

end

surf(XV,YV,ZV); shading interp; view(2); %colorbar; 
title('elevation')
figure();
surf(XV,YV,Z_hist); shading interp; view(2); %colorbar; 
title('history of the data')
%% transect data
if data_source(2)==1
    
    % url = 'http://micore.wldelft.nl/opendap/rijkswaterstaat/jarkus/transect.nc';
    url = 'http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/transect.nc';
    
    % readTransectDataNetcdf(url,'Ameland',[0 4953],2008);

    % nc_dump(url);

    year = nc_varget(url,'year');

    areacode = nc_varget(url,'areacode');
    ind = find(areacode==3);
    % indr = [min(ind) max(ind)];

    xj = nc_varget(url,'x',[min(ind)-1 0],[length(ind),1925]); 
    yj = nc_varget(url,'y',[min(ind)-1 0],[length(ind),1925]);
    zj = nc_varget(url,'height',[length(year)-1 min(ind)-1 0],[1,length(ind),1925]); 
     
    % xjall = nc_varget(url,'x',[min(areacode)-1 0],[length(areacode)-1,1925]); 
    % yjall = nc_varget(url,'y',[min(areacode)-1 0],[length(areacode)-1,1925]);
    % zjall = nc_varget(url,'height',[length(year)-1 min(areacode)-1 0],[1,length(areacode)-1,1925]);
    % figure(2); surf(xjall,yjall,zjall); shading interp; colorbar; view(2);
    
    [XJ,YJ]=ConvertCoordinates(xj,yj,'Amersfoort / RD New','xy','WGS 84','geo',CoordinateSystems,Operations); 
    [XJ2,YJ2]=ConvertCoordinates(XJ,YJ,'WGS 84','geo','WGS 84 / UTM zone 31N','xy',CoordinateSystems,Operations); 
    
    ZJ = zj;
    
    ind2 = find(abs(XJ(:,1)-5.8)==min(abs(XJ(:,1)-5.8)));
    figure(1); 
    subplot(211);
    surf(XJ2,YJ2,ZJ); shading interp; view(2); axis image; colorbar;
    subplot(212);
    plot(YJ2(ind2,:),ZJ(ind2,:),'r.');
end

%% Merge bathymetries from different sources

% interpolate jarkus data on grid vaklodingen
% ZJopV = interp2(XJ,YJ,ZJ,XV,YV);
% ZJopV = griddata(XJ,YJ,ZJ,XV,YV);

% make vectors xyz
% xyz = [XV(:) YV(:) ZV(:)];
% xyz{2} = [XJ2(:) YJ2(:) ZJ(:)];

%% create grid from XYZ-data

[X Y Z alfa VAR] = XBeach_GridOrientation(XV, YV,ZV,'dx',20,'dy',20,'manual','true');

XB = XBeach_selectgrid(X, Y, Z,'manual','false',CreateEmptyXBeachVar,VAR{:},'alfa', alfa,'posdwn', -1);

figure
pcolor(XB.Input.xInitial, XB.Input.yInitial, XB.Input.zInitial)
axis equal

%% write d3d and XBeach grid files...
if d3d==1
    ok=wlgrid('write','grid',XB.Input.xInitial,XB.Input.yInitial); 
elseif XBeach==1
    fi=fopen('bathy.dep','wt');
    for j=1:XB.settings.Grid.ny+1
        fprintf(fi,'%7.3f ',XB.Input.zInitial(:,j));
        fprintf(fi,'\n');
    end
    fclose(fi);

    fi=fopen('x.dep','wt');
    for j=1:XB.settings.Grid.ny+1
        fprintf(fi,'%7.3f ',XB.Input.xInitial(:,j));
        fprintf(fi,'\n');
    end
    fclose(fi);

    fi=fopen('y.dep','wt');
    for j=1:XB.settings.Grid.ny+1
        fprintf(fi,'%7.3f ',XB.Input.yInitial(:,j));
        fprintf(fi,'\n');
    end
    fclose(fi);

    fi=fopen('griddata.txt','wt');
    fprintf(fi,'nx      = %3i \n',XB.settings.Grid.nx);
    fprintf(fi,'ny      = %3i \n',XB.settings.Grid.ny);
    fprintf(fi,'dx      = %6.1f \n',XB.settings.Grid.dx);
    fprintf(fi,'dy      = %6.1f \n',XB.settings.Grid.dy);
    fprintf(fi,'xori    = %10.2f \n',XB.settings.Grid.xori);
    fprintf(fi,'yori    = %10.2f \n',XB.settings.Grid.yori);
    fprintf(fi,'alfa    = %10.2f \n',XB.settings.Grid.alfa);
    fprintf(fi,'depfile = bathy.dep \n');
    fprintf(fi,'vardx   = 1 \n');
    fprintf(fi,'xfile   = x.dep \n');
    fprintf(fi,'yfile   = y.dep \n');
    fprintf(fi,'posdwn  = -1 \n');
    fclose(fi);
end
    

    

