A few solutions come to mind that all take advantage of interpolation, depending on how 'ignore' is interpreted (I used NaN but you could use ==0 in place of isnan(z)):
% some fake data
[x,y] = meshgrid(-5:5,-5:5);
z = 50 - (x.^2 + y.^2);
nexttile
surf(x,y,z)
title('original')
% mark as nans, noting that the neighbor of every non-nan is nan:
z(1:2:end,:)=nan;
z(:,1:2:end)=nan;
nexttile
surf(x,y,z);
title('problem')
% soln 1: interpolate using nearest
badind=isnan(z);
f=scatteredInterpolant(x(~badind),y(~badind),z(~badind),'nearest');
zfixed=z;
zfixed(badind)=f(x(badind),y(badind));
nexttile
surf(x,y,zfixed)
title('soln 1')
% soln 2: interpolate using linear (smoother)
f=scatteredInterpolant(x(~badind),y(~badind),z(~badind),'linear');
zfixed=z;
zfixed(badind)=f(x(badind),y(badind));
nexttile
surf(x,y,zfixed)
title('soln 2')
% soln 3: interpolate to generate 4 values for every non-nan
nexttile
zi = interp2(z);
xi = linspace(x(1),x(end), width(zi));
yi = linspace(y(1),y(end), height(zi));
surf(xi,yi,interp2(z,'nearest'))
title('soln 3')
