function DM_startdigging

global data

if isempty(data.duneline)
    errordlg('No duneline detected. Please load/draw a duneline and try again.','No Duneline')
    return
end

data.znew = data.z;
if isnumeric(data.duneline);
    crest_line = data.duneline; %If crest_line was loaded from an external file
else
    crest_line = getPosition(data.duneline); %If crest_line was created using impoly
end
crest_line = interp1(crest_line(:,2),crest_line(:,1),data.y(:,1)); %Interp crest line 2 grid
Vcut_all = NaN(size(crest_line));
Vdune_all = Vcut_all;
yres = abs(nanmean(diff(data.y(:,1))));


%Perform preliminary duneline check to see if duneline is within elevation
%limits

I = find(~isnan(crest_line));
Zoriginal = interp2(data.x,data.y,data.z,crest_line(I),data.y(I,1));
if min(Zoriginal)<=data.startheight%||max(Zoriginal)>data.duneheight
    errordlg('Duneline is poorly located. Please change duneline position or dune dimensions and try again.','Poorly located duneline')
    return
end  

h = waitbar(0,'Calculating new grid....');
jj = 1;
for i = 1:length(crest_line)
    if isnan(crest_line(i))
        continue
    else
        waitbar(jj/length(I),h,'Calculating new grid....')
        jj = jj+1;
        [out,Vdune,Vcut] = make_artificial_dune([data.x(1,:)' data.z(i,:)'],crest_line(i),data.duneheight,data.dunewidth,data.startheight,data.duneslope_sea,data.duneslope_land,data.digoption);
        data.znew(i,:) = out(:,2)';
        Vdune_all(i) = Vdune;
        Vcut_all(i) = Vcut;
    end
end
close(h)

Vdune_total = nansum(Vdune_all)*yres;
Vcut_total = nansum(Vcut_all)*yres;
options.WindowStyle='replace';
options.Interpreter='tex';
msgbox({['Total dune volume = ' num2str(Vdune_total,'%4.0f') ' m^3'],['Total digging volume = ' num2str(Vcut_total,'%4.0f') ' m^3']...
    ,['Difference = ' num2str(Vdune_total-Vcut_total,'%4.0f') ' m^3']},'Dune statistics',options)
