% This script combines all the beach topography 
% data into one data cube with the same X and Y 
% coordinates and a selected spatial cell. 
% Each layer its a different field survey 

%% LOAD DATA  
clear
clc

disp('Indicate if they are data are from Cortadura (area=1) or from Camposoto (area=2)')
area=input('area=');

switch area
    case 1
        site = 'cort';
    case 2
        site = 'camp';
end



disp('Add information to an existing data-set (data=1) or create new data-set(data=2)')
data=input('data=');


dir = pwd;
f = filesep;
lchar = max(strfind(dir,f));
fw = dir(lchar+3:end);
switch data
    case 1   
        com_load = ['load ' dir(1:lchar) site '_alldata'];
        eval(com_load)
        [num,l]= size(who);
        fin = num-7;   
        name = [site num2str(fw)];
        comload =['load ' name '.mat'];
        eval(comload);
        comsave =['save(''' dir(1:lchar) site '_alldata'' ,''' site '*'')'];
        eval(comsave)
%        [sz]=size(ZALL,3);
    case 2
        fin =1;
        name = [site num2str(fw)];
        comload =['load ' name '.mat'];
        eval(comload);
        comsave =['save(''' dir(1:lchar) site '_alldata'',''' site num2str(fw) ''')'];
        eval(comsave)
end

    if area ==1
        load C:\MICORE\Cortadura\cortwall.mat
        corner = [max(wallcort(:,1)) max(wallcort(:,2)) 0];
        wallcort = [wallcort; corner];
    end

%if sz>fw
%    disp('THIS FILE IS BEEN GRIDED BEFORE')
%end

%% FIND THE  CORNERS OF THE GRID FROM THE FIRST 10 DATA-SETS
% if fin > 6
%     fin1 = 6;
% else
%   %  fin = fin+1;
% end

if fin<6
    for i=1:fin
        comXmin = ['Xmin(i) = min(' site num2str(i) '(:,1));'];
        eval(comXmin);
        comXmax = ['Xmax(i) = max(' site num2str(i) '(:,1));'];
        eval(comXmax);    
    
        comYmin = ['Ymin(i) = min(' site num2str(i) '(:,2));'];
        eval(comYmin);
    
        comYmax = ['Ymax(i) = max(' site num2str(i) '(:,2));'];
        eval(comYmax);
    end
    Xmin = floor(min(Xmin));
    Ymin = floor(min(Ymin));
    Xmax =  ceil(max(Xmax));
    Ymax =  ceil(max(Ymax));   
    
else 
        com_load1 = ['load ' dir(1:lchar) site '_allgrids'];
        eval(com_load1)
end

%% 
% Make the grid 

disp('Set the grid size (ds)')
ds=input('ds=');

if fin<6
    [x,y] = meshgrid([Xmin:ds:Xmax],[Ymin:ds:Ymax]);
else
end



if fin<6
    for i=1:fin
        if area ==1
            if i==1
                in = inpolygon(x,y,wallcort(:,1),wallcort(:,2));
            end
            comwall = [site num2str(i) '= [' site num2str(i) '; wallcort];'];
            eval(comwall)
        end

        comgrid = ['[z' num2str(i) ']=griddata(' site num2str(i) '(:,1),' site num2str(i) ...
           '(:,2),' site num2str(i) '(:,3),x,y);'];
        eval(comgrid)
        
         if area ==1
             comclean_cort = ['z' num2str(i) '(in)=NaN;'];
             eval(comclean_cort)
         end
        
        comclean= ['[z' num2str(i) ']=cleangrid(x,y,z' num2str(i) ',' site num2str(i) ');'];
        eval(comclean);
        
        comcomb = ['ZALL(:,:,' num2str(i) ')= z' num2str(i) ';'];
        eval(comcomb)
    end
else
    if area ==1
        comwall = [site num2str(fw) '= [' site num2str(fw) '; wallcort];'];
        eval(comwall)
    end
    
    comgrid = ['[z' num2str(fw) ']=griddata(' site num2str(fw) '(:,1),' site num2str(fw) ...
           '(:,2),' site num2str(fw) '(:,3),x,y);'];
    eval(comgrid)
        
    if area==1
        comclean_cort = ['z' num2str(fw) '(in)=NaN;'];
        eval(comclean_cort)
    end
        
    comclean= ['[z' num2str(fw) ']=cleangrid(x,y,z' num2str(fw) ',' site num2str(fw) ');'];
    eval(comclean);
        
    comcomb = ['ZALL(:,:,' fw ')= z' num2str(fw) ';'];
    eval(comcomb)   
end

if area ==1
    comsave1 =['save(''' dir(1:lchar) site '_allgrids'',''x'',''in'',''y'',''ZALL'')'];
    eval(comsave1)
else
     comsave1 =['save(''' dir(1:lchar) site '_allgrids'',''x'',''y'',''ZALL'')'];
    eval(comsave1)   
end
        

