% main file to turn NetCDF into a structure for topo analysis.

% needs Vlugtenburg TUD tools
% needs OE tools

% M.A. de Schipper 2009-2013

clear all
close all

%% Load
LoadSurveyData

%% Paths:  Rotate

%     x0=69141.4936;      %% coordinates northern Jarkus pole (exracted from the 5 Aug 2009 survey)
%     y0=447448.7293;

x0=67973; % Southern beach pole
y0= 446152;
theta0 =-131.65534968173412+90;  %% angle ref line (fence) in add_surveygrid!

for idate=1:length(pairs)
    X=Vsurvey(idate,1).paths.combinedpath(:,1);
    Y=Vsurvey(idate,1).paths.combinedpath(:,2);
    [Xr,Yr]=RotateTopo(X,Y,x0,y0,theta0);
    Vsurvey(idate,1).paths.combinedpath_rot=[ Xr Yr Vsurvey(idate,1).paths.combinedpath(:,3)];
end


%% Despike!
if 1
    load ([pathvlugtenburg 'Vsurveybadpoints.mat']);
    for idate=1:length(pairs)
        Vsurvey(idate,1).paths.combinedpath_rot(Vsurveybadpoints(idate,1).indexes,:)=[]; % remove bad spots
    end
end


% Maak transect data opnieuw

for idate=1:length(pairs)
    idate
    data.data=Vsurvey(idate,1).paths.combinedpath;
     data.data(Vsurveybadpoints(idate,1).indexes,:)=[];
 load surveylines
     % Now extract all surveylines
    for itransect=1:22
        transect= extract_line_nr_na_NetCDF(data,Surveylines{itransect});  %% do for all transects
        Vsurvey(idate,1).transects.elevation2(:,itransect)=transect.zi;
        
        
        disp(['   line ' num2str(itransect) ' extracted'])
    end
end





%% Make grids
dx=2;
dy=10;

% Unrotated
[Xi,Yi]=meshgrid(6.7E4:dy:6.95E4,4.46E5:dx:4.485E5);
for idate=1:length(pairs)
    Zi=griddata(Vsurvey(idate,1).paths.combinedpath(:,1),Vsurvey(idate,1).paths.combinedpath(:,2),Vsurvey(idate,1).paths.combinedpath(:,3),Xi,Yi) ;
    Vsurvey(idate,1).grids.Xi=Xi;
    Vsurvey(idate,1).grids.Yi=Yi;
    Vsurvey(idate,1).grids.Zi=Zi;
end

% Rotated
[Xi,Yi]=meshgrid(0:dy:1750,0:dx:1300);
for idate=1:length(pairs)
    Zi=griddata(Vsurvey(idate,1).paths.combinedpath_rot(:,1),Vsurvey(idate,1).paths.combinedpath_rot(:,2),Vsurvey(idate,1).paths.combinedpath_rot(:,3),Xi,Yi) ;
    Vsurvey(idate,1).grids.Xi_r=Xi;
    Vsurvey(idate,1).grids.Yi_r=Yi;
    Vsurvey(idate,1).grids.Zi_r=Zi;
end


%% Transects
for idate=1:length(pairs)
    [cross_shore_i,altitude_i] = combine_surveys_nr(survey_ind,idate);
    Vsurvey(idate,1).transects.cr_sh(:,:)=cross_shore_i;
    Vsurvey(idate,1).transects.elevation(:,:)=altitude_i;
end

save([pathvlugtenburg 'AllVbeachsurveys.mat'],'Vsurvey')
disp('saving done')

save([pathvlugtenburg 'datesandpairs.mat'],'dates','pairs')
disp('saving done')


%% Make Image

if 0
    idate=5;
    figure
    scatter3(Vsurvey( idate,1).paths.combinedpath(:,1),Vsurvey( idate,1).paths.combinedpath(:,2),Vsurvey( idate,1).paths.combinedpath(:,3),5,Vsurvey( idate,1).paths.combinedpath(:,3) )
    view(2)
    axis equal
end

%%
if 0
    figure
    for idate=1:length(pairs)
        %idate=5;
        pcolor(Vsurvey( idate,1).grids.Xi,Vsurvey( idate,1).grids.Yi,Vsurvey( idate,1).grids.Zi )
        shading interp
        caxis([-10 5])
        colormap(jet)
        colorbar
        pause(1)
        
        
    end
end


