close all; clear all;

sdir = pwd;

cd F:\subversion_tools\OpenEarthTools\matlab\

oetsettings

cd (sdir);

data_source = [1 1]

% addpath(genpath('F:\subversion_tools\McTools\mc_programs\DelftAlmighty\SuperTrans\'));
addpath(genpath('F:\subversion_tools\McTools\mc_programs\DelftDashBoard\general\SuperTrans\'));
load('CoordinateSystems.mat')
load('Operations.mat');

%% vaklodingen

if data_source(1)==1
    % TO DO: make generic function for easy use in future...
    % [X,Y,Z] = GetGridData(P)

    info = nc_info('warped2005.nc');
    % nc_dump(info.Filename);
    GT = str2num(info.Dataset(1).Attribute(6).Value); 

    GZ = [info.Dimension.Length];

    xv = (GT(1)+[1:GZ(1)]*GT(2));
    yv = (GT(4)+[1:GZ(2)]*GT(6));

    % use function GetPointsPolygon 
    indx = find(xv>5.62 & xv<5.99);
    indy = find(yv>53.40 & yv<53.53);

    Z = nc_varget(info.Filename,'Band1',[min(indy)-1 min(indx)-1],[length(indy),length(indx)]);

    [XV,YV] = meshgrid(xv(indx),yv(indy));
    
    [XV2,YV2]=ConvertCoordinates(XV,YV,'WGS 84','geo','WGS 84 / UTM zone 31N','xy',CoordinateSystems,Operations); 

    Z(Z==-32768) = NaN;
    ZV = Z/100;
    
    ind1 = find(abs(XV(1,:)-5.8)==min(abs(XV(1,:)-5.8)));
    figure(1); 
    subplot(211); 
    surf(XV2,YV2,ZV); shading interp; view(2); axis image; colorbar; hold on;
    subplot(212);
    plot(YV2(:,ind1),ZV(:,ind1),'b.'); hold on;
end

%% Jarkus raaien
if data_source(2)==1
    
    % url = 'http://micore.wldelft.nl/opendap/rijkswaterstaat/jarkus/transect.nc';
    url = 'http://dtvirt5.deltares.nl:8080/thredds/dodsC/opendap/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{1} = [XV2(:) YV2(:) ZV(:)];
xyz{2} = [XJ2(:) YJ2(:) ZJ(:)];

%% create grid from XYZ-data

[X Y Z alfa VAR] = XBeach_GridOrientation(xyz,'dx',100,'dy',100);

% [Xgr Ygr Zgr alfa VAR] = XBeach_GridOrientation(XV, YV, ZV,'manual','false','dx', 0.001,'dy', 0.001,'xori', min(Xarr),'yori', min(Yarr),'xend_y0', [max(Xarr) ],'x_yend', [0 max(Yarr)]);

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

% save d3dgrid
addpath F:\d3d\toolbox\
wlsettings('F:\d3d\toolbox\');

ok=wlgrid('write','grid',XB.Input.xInitial,XB.Input.yInitial); 



% close all;

%%
    

    

